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On the non-linear mechanics of hydrodynamic stability 


By J. T. STUART 
National Physical Laboratory, Teddington, Middlesex 


(Received 11 November 1957) 


SUMMARY 

In most work on the theory of stability of laminar flow, 
infinitesimal disturbances only have been considered, so that 
only the initial growth of the disturbance has been determined. 
It is the object of the present paper to extend the theory to larger 
amplitudes and to study the mechanics of disturbance growth with 
the inherent non-linearity of the hydrodynamical system taken 
into account. 

The Reynolds stress (where averages are taken with respect 
to some suitable space coordinate) is the fundamental consequence 
of the non-linearity, and its effects can be anticipated as follows. 
Initially a disturbance grows exponentially with time according 
to the linear theory, but eventually it reaches such a size that the 
transport of momentum by the finite fluctuations is appreciable 
and the associated mean stress (the Reynolds stress) then has an 
appreciable effect on the mean flow. ‘This distortion of the mean 
flow modifies the rate of transfer of energy from the mean flow 
to the disturbance and, since this energy transfer is the cause of 
the growth of the disturbance, there is a modification of the rate 
of growth of the latter. 

It is suggested that, in many cases, an equilibrium state may 
be possible in which the rate of transfer of energy from the 
(distorted) mean flow to the disturbance balances precisely the 
rate of viscous dissipation of the energy of disturbance. A theory 
based on certain assumptions about the energy flow is given to 
describe both the growth of the disturbance and the final 
equilibrium state, and application is made to the cases of 
Poiseuille flow between parallel planes and flow between rotating 
cylinders. ‘The distorted mean flow in the equilibrium state can 
be calculated and from this, in the latter case, the torque required 
to maintain the cylinders in motion. Good agreement is obtained 
with G. I. Taylor’s measurements of the torque for the case when 
the inner cylinder rotates and the outer cylinder is at rest. 


1. INTRODUCTION 
The immediate objective of the theory of hydrodynamic stability is to 
understand the mechanism of instability in laminar flow and to obtain a 
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criterion for its occurrence. A more fundamental objective is to understand 
how, and under what circumstances, turbulence may arise from laminar 
instability. ‘Ihe cpnnection between laminar instability and turbulence 
may not be a direct one, but under certain circumstances instability 
of laminar flow wil! be a necessary prelude to transition to turbulence. 
It is clear that the stability problem in its general form must be 
considered to be non-linear, because the equations of motion are non- 
linear. 

The mathematical problem of hydrodynamic stability can be formulated 
by taking the given steady-state solution of the equations of motion and 
superimposing a disturbance of a suitable kind; this results in a set of 
(non-linear) ‘ disturbance’ equations, which govern the behaviour of the 
disturbance. If the solution of the equations shows that any disturbance 
ultimately decays to zero, the flow is said to be stable; whereas if the 
disturbance can be permanently different from zero, the flow is unstable. 
It does not always happen that instability leads to turbulent motion, because 
another (possibly more complex) form of laminar motion may be the result. 
Indeed, it will be shown that this is often the case. 

Naturally, the solution of the disturbance differential equations is 
simplified considerably by linearization for small disturbances, and for 
a description of theories based on this assumption the reader is referred 
to the book by Lin (1955). On the basis of linear theory it is possible to 
consider disturbances which contain an exponential time factor of the form 
exp(kt), t being the time. ‘The boundary conditions on the disturbance 
equations require the vanishing at the boundaries of quantities like the 
disturbance velocity components. Consequently, the boundary conditions 
are homogeneous, and there is an eigenvalue problem for the determination 
of the quantity &. In this (linear) case, stability or instability is defined 
as follows: if it is possible for k to have a positive real part, the flow is 
unstable; otherwise the flow is stable. 

‘The prediction by linear theory of a disturbance which increases 
exponentially with time is a feature which has occasionally given rise to 
the suggestion that turbulence would necessarily ensue from the growth 
of the disturbance to large amplitudes; however, examples are known 
in which this is not the case. On the other hand, it has been argued that 
the non-linear terms will stabilize completely a flow which is unstable 
according to linear theory, but such arguments can generally be refuted 
(Stuart 1956a). In this paper, certain features of the role played by the 
non-linear terms of the equations of motion are discussed, and this leads 
to a clarification of the connection between linear and non-linear instability 
theories. Attention is restricted to flows which have constant local 
Reynolds number. 

In cases of instability of fluid flow, the disturbance is usually periodic 


in at least one direction of space. ‘Thus it is convenient to take averages 
with respect to one of the spatial dimensions, and to separate the flow into 
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a mean part and a disturbance part, where the latter has zero mean*. It is 
clear that the two parts of the flow are interdependent through the action 
of the Reynolds stress (arising from the disturbance) on the mean flow. 
(For a discussion of Reynolds stresses, the reader is referred to the books 
by Goldstein (1938) and ‘Townsend (1956).) On the basis of linear theory, 
the disturbance is assumed to be so small that the effect of the Reynolds 
stress on the mean motion can be neglected, in which case the mean flow 
is the original laminar flow. However, in a non-linear theory the inter- 
dependence of the mean and disturbance parts of the flow must be taken 
into account. Let us now consider a flow whose local Reynolds number 
does not vary, as in the case of flow between parallel planes or concentric 
cylinders, and let the flow be perturbed by a small disturbance. Initially 
the disturbance amplifies exponentially with time according to linear theory, 
but eventually it reaches such a size that the mean transport of momentum 
by the finite fluctuations is appreciable, and then the associated mean stress 
(the Reynolds stress) has an appreciable effect on the mean flow. This 
distortion of the mean flow clearly modifies the rate of transfer of energy 
from the mean flow to the disturbance and, since this energy transfer is 
the cause of the growth of the disturbance, there is a modification of the 
rate of growth of the latter. ‘These processes, in which the disturbance 
distorts the mean flow and the distortion of the mean flow modifies the rate 
of growth of the disturbance, occur simultaneously. 

It is natural to enquire if an equilibrium state is possible, in which the 
rate of transfer of energy from the mean flow to the disturbance balances 
precisely the rate of viscous dissipation of energy of the disturbance. In 
such an equilibrium state, the disturbance will have a definite finite amplitude 
and the mean flow will be distorted from its original laminar form. 
Experimental evidence of an equilibrium state of this kind is afforded by 
G. I. Taylor’s observations on the instability of flow between rotating 
cylinders, where the instability takes the form of cellular, toroidal vortices 
spaced regularly along the axes of the cylinders. ‘Taylor (1923, p. 342) 
observed that ‘‘ A moderate increase in the speed of the apparatus merely 
increased the vigour of the circulation in the vortices without altering 
appreciably their spacing or position ’’, and suggested that *‘ ‘he experiments 

. indicate that the effect of the second-order [non-linear] terms is to prevent 

* It should be mentioned that, an another formulation of the problem, the ‘ dis- 
turbance ’ is defined to be the whole of the deviation from the original laminar flow. 
However, with this definition the ‘ disturbance’ must contain a ‘ mean’ part, a fact 
which was overlooked by several authors in their studies of non-linear instability 
theory (Stuart 1956 a). In the opinion of the writer, this approach does not yield 
such a clear understanding of the physical processes involved in instability as does the 


‘ 


approach based on the concept of a mean flow and a disturbance flow (with zero mean) 
interacting through the action of a Reynolds stress. Consequently, in this paper the 
flow will always be separated into a mean part and a disturbance part (with zero mean). 
The term ‘ disturbance’ will not be used to denote the deviation from the original 
laminar flow. 
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the vortices from increasing indefinitely in activity’’. “The combination of 
the mean flow with a disturbance of definite amplitude may be referred to 
as an equilibrium flow. Examples relating to Poiseuille flow between parallel 
planes and to flow between rotating cylinders are described in §2 and §3, 
and for the case of flow between cylinders good agreement with experiment 
is obtained. For the case of Poiseuille flow, however, there is no experimental 
evidence that an equilibrium flow of the kind described above does occur. 
There is the possibility in this case, as in other cases, that a flow of this kind 
does not occur, but rather that there is a continual generation of harmonics 
of the basic disturbance and of other disturbances. 

In the discussion above, attention has been paid to equilibrium flows 
which may develop when the original laminar flow is unstable according to 
linear theory. Some flows, however, such as Couette flow between parallel 
planes and Poiseuille flow in a circular pipe, are completely stable against 
infinitesimal disturbances. Even so, turbulence can occur at sufficiently 
high Reynolds numbers. Furthermore, turbulence occurs in some flows 
(for example, in Poiseuille flow between parallel planes) at a lower Reynolds 
number than the critical according to linear theory, that is, it occurs when the 
laminar flow is stable with respect to infinitesimal disturbances. A suggestion 
which may lead to an explanation of such phenomena is that the appropriate 
laminar flow may be unstable with respect to finite disturbances. When a 
disturbance of suitable magnitude is present, the mean flow may be 
distorted to such a form that the rate of transfer of energy to the disturbance 
can balance exactly its rate of dissipation by viscosity. On the other hand, 
a finite disturbance which is small enough will presumably decay to zero 
amplitude, either because the rate of energy transfer to the disturbance is 
insufficient to balance the rate of viscous dissipation of kinetic energy, or 
because the energy transfer is actually from the disturbance to the mean 
flow. As an example of instability for finite disturbances, Meksyn & Stuart 
(1951) considered the case of Poiseuille flow between parallel planes and 
showed that the critical Reynolds number drops as the amplitude of the 
disturbance rises. ‘The reader is referred also to Lin (1955) and Stuart 
(1956b) for a discussion of these ideas. A discussion of the formulation of 
the instability problem for plane Couette flow with finite disturbances was 
given by Noether (1921). 

It will be convenient to refer to non-linear disturbances as existing under 
supercritical conditions if the Reynolds number is above the value which is 
critical for linearized instability, and as existing under subcritical conditions 
if the Reynolds number is such that the flow is stable with respect to 
infinitesimal disturbances. A non-linear disturbance may clearly arise 
spontaneously under supercritical conditions as a result of continued 
amplification of a secondary disturbance; the sequence of events which 
would lead to a non-linear disturbance under subcritical conditions is less 


evident and will not be discussed here. 
If equilibrium flow is established consisting of a mean flow with a 
steady-amplitude finite disturbance, it does not follow that it is a stable flow 
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at all Reynolds numbers above the critical value for which it first occurs. 
It is likely, however, to be stable for a certain range of Reynolds number 
above the critical. Moreover, it is clear that it must become unstable at 
some Reynolds number if turbulence is to occur at higher Reynolds numbers. 
For example, in the case of flow between rotating cylinders mentioned 
above, ‘Taylor observed that ‘‘a large increase [in the angular speed] 
caused the symmetrical motion [of the cellular vortices] to break down 
into some kind of turbulent motion, which it was impossible to follow by 
eye’’. Another interesting case is that of Poiseuille flow between parallel 
planes, where the critical Reynolds number according to linear theory is 
a Reynolds number for which the flow is normally turbulent. Consequently, 
equilibrium flows under supercritical conditions are, in this case, almost 
certainly unstable. Similar considerations apply to equilibrium flows under 
subcritical conditions, at least for the higher Reynolds numbers for which 
they are valid. 

‘The development of a non-linear instability theory for boundary-layer 
flows, where the flow and local Reynolds number change in the stream 
direction, presents additional difficulties. Whereas in the linear instability 
theory it is permissible to regard a boundary-layer flow as nearly parallel 
and to neglect boundary-layer growth, it does not seem obvious that such 
an approximation is permissible in a non-linear theory. Because the local 
Reynolds number increases in the downstream direction, any disturbance 
is convected into regions of higher Reynolds number and the effect of this 
continuous change of Reynolds number would have to be taken into account. 
Consequently, the non-linear theory which will be described in this paper 
for flows with constant local Reynolds number does not necessarily apply 
quantitatively to the case of the boundary layer. ‘he main features of the 
theory are of wide applicability, but there may be additional factors 
influencing them because of the growth of the thickness of the boundary 
layer. 

An interesting suggestion concerning the development of turbulence 
from the growth of small disturbances has been made by Landau (1944). 
As noticed above, the occurrence of instability in a flow may lead to the 
replacement of the original laminar flow by a new laminar flow, which consists 
of a mean flow with a superimposed finite disturbance. ‘This flow may 
be expected to persist for a certain range of Reynolds number above the 
critical value and then to become unstable at some Reynolds number against 
a new (second) type of disturbance. A new equilibrium flow, consisting 
of a mean flow with two superimposed modes of disturbance, is then 
conceivable for a range of Reynolds number above the second critical value. 
As the Reynolds number is raised still further, additional modes of 
disturbance may appear successively until, at sufficiently large Reynolds 
number, the flow is so highly disturbed as to be considered turbulent. In 
the case of flow between rotating cylinders, experiments show that the 
development of turbulence takes place fairly slowly as the Reynolds number 
is raised; this would correspond to the first two at least of the successive 
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critical Reynolds numbers of Landau’s theory being fairly widely spaced. 
On the other hand, there are cases of flow in which turbulence develops 
rather suddenly as the Reynolds number is raised, and in these cases one 
might infer that the critical Reynolds numbers are close together. 

In the discussion given above, attention has been entirely directed to 
instability in fluids which are in motion. ‘There is, however, another 
important type of instability, namely that which occurs when a horizontal 
layer of fluid is heated from below, causing an unstable density gradient 
in the fluid. Instability takes the form of convection in cells of polygonal 
planform, and occurs when a certain parameter (the Rayleigh number) 
is above a critical value. ‘The linear theory of instability is well understood 
(see, for example, Pellew & Southwell 1940) and is known to be analogous, 
under certain conditions, to the theory of instability of flow between 
rotating cylinders. ‘The fundamental physical process in the instability 
is the conversion of potential energy associated with the gravitational field 
into kinetic energy of the convective disturbance motion. When the 
disturbance has a finite amplitude, the mean heat transport by the 
convective motion causes a modification of the mean temperature distribu- 
tion, where averages are taken over horizontal planes; therefore, because 
of the corresponding modification of the mean density distribution, the rate 
of transfer of energy from the gravitational field to the disturbance is 
modified. It appears that a steady (equilibrium) state is possible, in which 
the rate of transfer of energy into the disturbance balances precisely the rate 
at which energy is dissipated and diffused. ‘Vhe non-linear theory of thermal 
convection has been studied by Sorokin (1954) and by Malkus & Veronis 
(1958), both by methods related to those of the present paper and by 
perturbation series expansions. 

Some of the results obtained in this paper are similar to results obtained 
from the turbulence models due to Burgers (1948) and Hopf (1956). For 
a discussion of non-linear instability effects in the Burgers model, the reader 
is referred to a paper by Stuart (1956b). 


2. DiIsTURBANCES UNDER SUPERCRITICAL CONDITIONS IN POISEUILLE FLOW 
BETWEEN PARALLEL PLANES 

In this section a simple treatment of the non-linear problem for 

disturbances under supercritical conditions is described, and it is applied 

in this section and in $3 to the cases of Poiseuille flow between parallel 

planes and of flow between rotating cylinders. The method is based 

essentially on a balance of energy, taking into account the mean-flow 


distortion due to the Reynolds stress, and it is the ‘integral’ properties 
of the flow which are treated rather than the spatial details. Consequently, 
it is clear that the method has an obvious application to cases of instability, 
such as that of flow between rotating cylinders, where it is the overall 
properties of the flow (as opposed to the details) which are impeztant. 
On the other hand, in cases of instability which are governed by the 
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Orr—Sommerfeld equation (for example, the case of plane Poiseuille flow), 
the instability depends very much upon precise details of the mean flow, 
especially in the vicinity of the critical layer. Consequently, the method 
is less valid for such cases, and the application to plane Poiseuille flow has 
less validity than the case of flow between rotating cylinders. Moreover, there 
is less experimental evidence of quantitative value in the former case than 
in the latter. In spite of these facts, the application to plane Poiseuille flow 
is of considerable illustrative value, and is therefore presented in this section 
before a treatment in $3 of the case of flow between rotating cylinders, for 
which a comparison with experiment will be made. 

Consider a two-dimensional flow between parallel planes, and let x, 
denote the coordinate parallel to the planes and x, the coordinate normal 
to them, with u,, u, denoting the corresponding components of velocity, 
Let P denote pressure, p density, v kinematic viscosity and ¢ the time. 
Now let us introduce a reference length L, equal to half the distance between 
the planes, and a reference speed U,, the maximum speed in laminar flow. 


Thus we define 
x, = Lx, x, = Lz, u,y= Uy yu, ua= Uw, P=ppU2Z, t=Lt/U,. (2.1) 


Then the Navier-Stokes equations and equation of continuity are 


Cu Cu ou C | 
——— i Se Ms a se V7u, 
ct Ox oP Ox R 
ow ou Cw 0 i 7 
— +u— +w—=-~+ Vm, } (2.2) 
ct Ox oP oh R 
Cu out 
— . a 0, 
CX Cog J 


where V? = 0?/dx? + 62/03", and R = U, L/v denotes the Reynolds number 
for the problem. 
The solution for undisturbed laminar Poiseuille motion is 


a= 1—22 w= 0 p = —2x/R+const., (2.3) 


@ 39 


where a bar above a quantity denotes a mean value with respect to the 
distance x. This is the basic laminar motion whose non-linear instability 
characteristics we wish to examine. 
Now let 
u= i(z,t)+u’' (x, 2, t), w= w'(x,2,t), (2.4) 


where uw’ and w’ denote components of a finite disturbance (with zero mean). 
Since the disturbance is finite, @ is no longer given by (2.3) but is distorted 
by the Reynolds stress. Furthermore, if the disturbance is growing or 
decaying in amplitude, 7 depends on time. We suppose that uw’ and w’ 
are harmonic in x with wavelength 27/x% and are given by Fourier series. 
If we substitute (2.4) into (2.2), separate out the mean and disturbance 
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parts of the motion and then integrate the resulting disturbance equations 
to obtain the energy balance relation, we have (without approximation) 
( { dae a a: ee Cu 1 rr few’ cu’ \2 
— || 3(u’2+w’*) dxdz || (—u’w’)— dxdz— — || (— — —] dxdz, 
CT ; ‘ COs R BE a \ , 

I,-1,/R, (2.5) 


where the integrands are evaluated over a volume bounded by the planes 





ui 


and by one wavelength. ‘The mean velocity @ occurring in (2.5) is given 
by a a A A = 
’ j 1 ca 
— T — = = _— T 7 a apr . ( 2 . 6 ) 
OT oz ox R oz 


For a discussion of the derivation of the above equations, and for the 


non-linear disturbance differential equations for w’ and w’, the reader is 
referred to Stuart (1956a,b). Equations (2.5) and (2.6) will be used 
here as an approximate basis for the solution of the non-linear problem 
of the growth of the disturbance (u’, w’). 

We note that in equation (2.5) the term on the left-hand side gives the 
rate of growth of the disturbance energy within the volume considered. 
On the right-hand side, the term J, is the integral of the product of the 
Reynolds stress and the mean velocity gradient, and represents the rate 
of energy transfer from the mean flow to the disturbance; the term J,/R 
is always positive and represents the rate of viscous dissipation of energy 
of the disturbance. If /, is positive and greater than J,/R, the disturbance 
energy is growing and the disturbance is increasing in amplitude. 
Equation (2.6) shows how the distribution of mean motion is affected by 
the viscous stress, pressure gradient, and Reynolds stress due to the 
disturbance. An equilibrium flow is possible if @ can be so distorted by 
the Reynolds stress that J, = /,/R for a given Reynolds number. 

It is worth considering at this point the iene: conditions on the 
mean motion. For the velocity they are @ = 0 at s + 1. We shall also 
ado pt the condition that the mean pressure gri adient shi lr remain unchanged; 
this is clearly physically realistic, since the pressure gradient is externally 
applied*. As a consequence, the mean skin friction on the walls will be 
the same in the coreg flow as in the laminar flow. However, the mean 
velocity near the channel centre must drop when a disturbance is present, 
because energy is required to maintain the disturbance. ‘The unchanged 
pressure gradient is given by equation (2.3) and the mean motion equation 
thus takes the form 

Cu C ——— 2 1 Cn al 
— +S = w ne peer Bs ge (2.7) 
ot dz m Kids" 

‘ 


In a state of equilibrium, @ is independent of time and equation (2.7) is 





easily integrated to yield 


u(z)= 1-—3s?+R | ww’ dz. (2.8) 


* Other boundary conditions, such as that of constancy of mass flux, yield qualita- 
tively similar results for the amplitude of oscillation. 
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Since u’w’ is an odd function of z, 7=0 at s = +1. Moreover, since 
I, is positive, u’w’ must be of opposite sign to éa/dz in at least a dominant 
range of 3; therefore, it must be dominantly positive when s is positive. 
Thus @ is less than unity at z = 0 and, in general, a(z) is everywhere less 
than its value for laminar flow. A corollary of the condition on the mean 
pressure gradient is that the Reynolds number is based on a velocity Uy 
which is the maximum for laminar flow with the same mean pressure gradient. 

Suppose we now define a stream function for the disturbed flow in the 


form 


~ 
Zialar—c, 7) 


wb = d,(2,t)+ $, (2,7 )e* sii al, (3,7 )e Oe eer (3, Te” 
+¢$.(z, BE cg CE 
where the symbol ~ represents a complex conjugate. uw represents a 
mean flow together with a periodic disturbance consisting of the first 
harmonic with wavelength 27/x, and higher harmonic components having 
wave-numbers nz (m integral) but the same (real) wave velocity c,, which 
is assumed to be independent of time. ‘The amplification or damping of 
the disturbance, and the consequent changes in the mean velocity (given by 
i = 0¢,/ ez), are accounted for by the dependence of all the ¢-functions on r. 
If we substitute (2.9) into (2.2), utilize 
u = cd/oez, w= — od/ex, (2.10) 
eliminate the pressure and separate out the Fourier components of the 
disturbance, we obtain a set of equations for the functions #@, 4,, ds, and 
so on. ‘The equations are similar to (3.2), (3.3) and (3.4) of the paper by 
Stuart (1956b), except that allowance must be made for the dependence 
of the functions on 7. Because of the non-linearity of the system, all of 
the functions are mutually dependent and an infinite set of equations has 
to be solved to obtain an exact solution. In particular, the mean-motion 
equation is (2.6) with 

uw = inh dy = id, T 2(¢, dy — $, ds) scot, (2.1 1) 

where primes denote differentials with respect to 3. 

In the linear theory of instability, the Reynolds stress is neglected, 
so that the mean velocity is the laminar Poiseuille flow. Furthermore, 
the stream function (2.9) contains only the functions ¢) and ¢,, the real 
part of which determines the actual stream function. ‘The problem then 
is the solution of the Orr-Sommerfeld equation (see Lin 1955), namely 

(a —c,—1¢;)(b| — #°6,)— a", = — (t/aR)(P¥ - 2476, +a4p,), (2.12) 
where @ = 1—2 and the function 4, is proportional to exp(zc¢;7). ‘The 
non-linear problem posed here is to follow the growth of a disturbance 
at a given Reynolds number and, in particular, to find out if an 
equilibrium flow is possible at that Reynolds number. ‘To this end, 
the stream function (2.9) has been chosen to allow for changes in the basic 
disturbance, partly through the generation of harmonics, due to the 
non-linearity. 
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A simple approximate method of solution 

The non-linear effects of a disturbance of finite amplitude appear to 
be of two kinds. On the one hand there is, through the Reynolds stress, 
a modification of the energy transfer between the mean flow and the first 
harmonic component, ¢,, of the disturbance. On the other hand, higher 
harmonics (d,, etc.) of the disturbance are generated and there is an 
interchange of energy between them and both the mean flow and the first 
harmonic component. ‘The assumption which will be made here is that the 
dominent non-linear interaction is that between the mean flow and the first 
harmonic component of the disturbance. ‘To put it in another way, the 
effect of the Reynolds stress on the mean flow is taken into account, thus 
causing modifications of the energy transfer between the mean flow and the 
tirst-harmonic disturbance, but the generation of higher harmonics is 
ignored. A similar assumption was made by Meksyn & Stuart (1951) 
for disturbances under subcritical conditions. In a formulation of this 
kind there are two differential equations to determine &@ and ¢,, namely, (2.7) 
and (2.12) with c; replaced by x~!¢/ez; in equation (2.11), d, = ¢, =... =0. 

In the present case of disturbances under supercritical conditions, we 
shall make the additional assumption that the disturbance (u’, w’) is 
similar in shape to a solution given by linear theory, but that the solution 
is multiplied by an amplitude factor, a(7), which is a function of time. 
Then the approximate problem is to satisfy the disturbance energy-balance 
equation (2.5) and the equation of mean motion (2.7), in which case the 
amplitude is determined. If the amplitude of the disturbance in an 





equilibrium flow is required, @ is given in terms of uw’ by (2.8) and can 
be substituted into (2.5). With the ‘shape’ assumption, (2.5) then gives 
a relation between a? and a’, so that a? is determined in terms of integrals 
involving functions given by linear theory. ‘This is the most important result. 

An equation illustrating the growth of the disturbance may be obtained 
by noting that the term ca#/cr in (2.7) is negligible at small times on the 
basis of linear theory, and also at large times because it is to be expected 
that a steady mean flow will be approached. Ii, on this basis, the term 
cu/er is ignored at all times, (2.5) and (2.8) yield the equation 

ae 
oe y, 2a" —y.%7Ra! a. (2.13) 
= R 

where a is the amplitude of the disturbance and 
1 l 
| S'f' |? +0716 2b dz, Vo 41 s(d db. d, db.) dz, 


(2.14) 
4| (¢,¢6,-4,46:')? ds, yy=2| {6"2+22? 46'?24+2!467) dz, 
0] 0 


the suffixes r and 7 denoting real and imaginary parts. The function ¢(z) 
y 


is the amplitude distribution of a disturbance stream function according 
to linear theory. Equation (2.14) is of the form 
dy 


— =f,v— Bov—Bs3y 
Pee ee 
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where y is the square of the amplitude of the disturbance. Equation (2.15) 
was given earlier by Landau (1944), although he does not state how, and 
under what assumptions, it was derived. 

The terms in (2.15) and the corresponding ones of (2.13) have the 
following physical meanings. The term dy/d7 represents the rate of change 
of the kinetic energy of the disturbance, the term §,y the rate of energy 
transfer from the mean flow to the disturbance, the term f,¥ the rate of 
viscous dissipation of the disturbance, and the non-linear term f, y* the 
restriction of the rate of energy transfer to the disturbance by the 
Reynolds-stress distortion of the mean flow. For very small disturbances, 
the non-linear term is negligible, and if 8, > 8, the disturbance amplifies 
like C exp[(8,—2)7], where C is a constant. ‘This solution corresponds 
to linear stability. On the other hand, if the Reynolds stress is included 
the general solution of equation (2.15) is 





Ce i 
y= - , (2.16) 
De ma = 
1+ { ——— }Ce's— Are 
Q QO 
Py— P2 
When 7 — ~, y~>Cexp[(8,—83)7], which is the linearized solution 
given above. Furthermore, when t— + ©, y+>(8,—[.)/x83, whatever 


the value of C. This solution suggests that, whatever the initial size of 
the disturbance, it always builds up to the same limiting amplitude. This 
conclusion is physically reasonable, although it should be remembered 
that the exact form of (2.16) is based on the assumption that ea#/er may 
be neglected at all times. 

Reverting now to (2.13), we note that the critical Reynolds number is 


R. = 4/472 


} 


the amplification factor of ¢? on linear theory is 


Y1 
and the square of the equilibrium amplitude is 
2 you—-y,/R  y(R-R,) 


a= = = = 


P3 ysaR yg 2R? 





Thus a, is proportional to the square root of the difference between the 
actual Reynolds number and the critical Reynolds number. In order to 
calculate specific values of the amplitude at a given Reynolds number, it is 
necessary to approximate to the function ¢. A fairly obvious choice would 
be the function appropriate to the minimum critical Reynolds number. No 
calculation of this function has been made, although Thomas (1953) has 
calculated and tabulated ¢ forx = 1, R = 104,¢ = c,+ic; = 0-2375 +0-00371. 
The function is illustrated in figure 1, # denoting the real part and / the 
imaginary part. Numerical integration of (2.14) then gave 


y, = 2:05146, yy. = 0-040192, , = 0-002308, y, = 247-62. (2.17) 


72 73 
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Using these values, we find that R= 6150, which isin reasonable agreement 
with Thomas’s critical value, namely, R, = 5780 at « = 1-02. ‘Thus, even 
though the function 4 is not the correct stream function at the critical 
Reynolds number, the energy-balance relation yields a fairly good approxi- 
mation to the critical Reynolds number. It is suggested, therefore, that 
the numbers (2.17) may be applicable over a wide range of R, and possibly 
also for small variations of z. Another characteristic of the linear instability 
theory which can be calculated from (2-17) is the amplification constant ¢, 
at « = 1, R= 10', and this is found to be c, = 0-00376. ‘The closeness 
of this to ‘Thomas’s value, 0-0037, serves as a check on the accuracy of the 


numerical work. 
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Figure 1. Function ¢ for Poiseuille flow at R = 104, a 1, c = 0:2375 +0-0037:. 


For the non-linear properties, we consider a wave-number, « = 1, 
since this is close to the value at the critical Reynolds number. The 
equilibrium amplitude number, a,, rises from zero at the critical Reynolds 
number (6150) through 0-00256 at R = 10! to a maximum value of 0-00266 
at R= 12300 (=2R.). To obtain the actual (dimensionless) stream 
function at a given value of R, the curves in figure 1 should be multiplied 
by the appropriate value of a. The distribution of Reynolds stress is shown 
in figure 2 (obtained from Thomas’s calculations), from which the actual 
(dimensionless) Reynolds stress can be obtained by multiplying by a’. 

Using the calculations described above, the mean velocity profile of 
the equilibrium flow can be calculated from equation (2.8); for R = 104, 
it has the form shown in figure 3. The distorted profile is everywhere less 
than the laminar profile because of the energy transfer to the disturbance, 
but has the same gradient at the wall. It can readily be shown that the work 
done by the pressure gradient to maintain the motion precisely balances the 
total viscous dissipation of energy in the equilibrium flow. The root-mean- 


squares of the disturbance velocity components, (u’?)!? and (w’?)!?, are 
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Figure 2. Reynolds stress function for Poiseuille flow at R = 10%. 
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Figure 3. Comparison of laminar and distorted mean flows. 
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shown in figure 4 for « = 1, R= 104. The magnitudes are similar to those 
applicable in turbulent motion at the same Reynolds number. A further 
point worthy of note is that, at a Reynolds number of 10*, the maximum 
value of the Reynolds stress is of the same order as the maximum value of 
the viscous stress. 

As mentioned in the Introduction, the equilibrium flow under super- 
critical conditions in Poiseuille flow between parallel planes is almost 
certainly unstable, because turbulence normally exists at those Reynolds 
numbers for which the supercritical flow exists. However, it is of interest 
that the magnitude (see figure 4) of disturbance which can be sustained 
is similar to the magnitude appropriate to turbulent flow at the same 
Reynolds number (Reichardt 1938; Laufer 1950). 
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Figure 4. Root-mean-square velocity fluctuations at R = 101. 


In connection with the establishment of the equilibrium flow, it has been 
suggested here that an infinitesimal disturbance at a given Reynolds number 
can amplify into an equilibrium flow. Experimentally, however, it is the 
equilibrium flow itself which is important, rather than the way in which 
itisattained. In practice the equilibrium flow may arise as follows. Suppose 
the flow can be kept free of disturbances, so that disturbances and turbulence 


do not occur under subcritical conditions. ‘Then, as the Reynolds number 
is raised, a weak disturbance will appear at the critical Reynolds number. 
The disturbance will have a small but finite amplitude. Then, as the 
Reynolds number is increased, the flow could (if it were not unstable) 
retain the same wavelength as the disturbance which first occurred, but 
with an amplitude a ~ R-1(R—R,)'*. This idea is discussed again in §3 
in connection with the flow between rotating cylinders. 





i 
; 
f 
: 





w 


On the non-linear mechanics of hydrodynamic stability 1 


3. DISTURBANCES UNDER SUPERCRITICAL CONDITIONS IN FLOW BETWEEN 
ROTATING CYLINDERS 
‘The second application of the theory will be to the flow between two 
concentric rotating cylinders, with 7, 6, s as the cylindrical coordinates 
and u, v, w the corresponding velocity components. It is assumed that 
the flow has rotational symmetry and is therefore independent of @. ‘Then 
the Navier-Stokes and continuity equations are 











— = ao = > 
ct or C38 p oT a 
Cz Ov ov Ut 1 
av t+u— +Ww— — ( V2 = =) 9 
( or Cs 7 Ys 
: ; : ~ (3.1) 
ou ou ou C ’ 
—— += +0 —— : Le 2 + vV2w, 
ct c} robs p cz 
1 oO ou 
— —(ru) 0, 
) Os 
where - oC 1 0 oC 
si Ge! = aS aes (3.2) 
oF Yr oj oZ2* 


It is known from the work of ‘Taylor (1923) and others that, when the 
rotational speed of the inner cylinder (or, to be more specific, a parameter 
sometimes called the ‘Taylor number) is above a critical value, the steady 
laminar flow is unstable. ‘The disturbance which appears takes the form 
of cellular, toroidal vortices spaced regularly along the axis of the cylinders. 
The linear instability theory for this flow is well known (‘Taylor 1923) 
and our purpose here is to study the non-linear theory, primarily in order 
to obtain the amplitude of the equilibrium flow. 

Since the flow which results from the instability is periodic with respect 
to s, it is convenient to take averages with respect to x. In order to allow 
for the distortion of the disturbance by the non-linearity, we write 
u=u' = u,(r, the + u(r, te? +...+4,(r, tle" +0,(7, the +... (3.3) 
together with exactly similar series for v’ = v—d@(r,t) and w’ = w, where 
the symbol ~ denotes a complex conjugate and a bar above a quantity 
denotes a mean value. If the above series are substituted into the equations 
of motion (3.1) and the Fourier components are separated, there results 
a set of equations for all the velocity functions involved in (3.3). 


For the mean motion, we find 


ra 1 aoe _ 1 0 > \ 
~ < (wt)— = (07+ 0) = —-Z, (3.4) 

Vr ¢ } po} 
ov ] oO yr o2 1 O 1 = 5 = 
aa *? te ar wn 


similar equations have been derived and discussed by ‘Townsend (1956) 
in connection with turbulent flow. Equation (3.4) gives the radial pressure 
gradient required to balance the centrifugal force and Reynolds stress, 
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gives the mean rotational velocity @ as a function of the Reynolds 
e dependence of @ on ¢ is retained because, if the disturbance 


while (3. 


>) 
stress. Th 
is growing or decaying, the mean motion undergoes distortion in order 
to maintain the energy balance. A suitable boundary condition to apply 
to the mean flow is that just enough external power is supplied to maintain 
the angular velocities of the cylinders at constant values, even though the 
skin friction is changing. If r, and r, are the radii of the inner and outer 
cylinders and @, and , their angular velocities, the boundary conditions 
are 
6=7,Q, atr=7n,; G=7,0), at r= Fz. (3.6) 
In a state of equilibrium, 0@ ct is zero and (3.5) can be integrated to give 


3 rir uy = 
@ = Ar+4 ae | ile: (3.7) 
. sie 7 
where __ 
(),(1 — mr2/r?) 1 Te y'y’ 
A= 2! +5 — | — ob, 3.8 
1 —r3/r? v(rz/rs—1) J, 4 i 
r?O.(1—m) Te uv’ 
B D ens! aie’ peerenweesl lr, 3.9 
F 1 —73/r5 Al —rVP8y J r ” (9.7) 


and m = Q,/Q,. With zero Reynolds stress, these relations become those 
appropriate to laminar Couette flow between rotating cylinders. 

In addition to the mean motion equations, it is possible to write down 
an infinite set of differential equations for the harmonic components of 
the disturbance. ‘These components are all mutually dependent because 
of the non-linearity of the system. ‘These equations will not be derived 
or considered in detail here, because it is proposed to use an approximation 
similar to that used in §2. ‘he assumption is that the dominant non-linear 
interaction is between the mean flow (@) and the first harmonic component 
of the disturbance (u,,7,,#,), where the latter is assumed to be similar 
in shape to that given by linear theory but multiplied by an amplitude 
factor. ‘he mean flow is distorted by the Reynolds stress of the disturbance 
and the consequent alteration of the energy transfer between mean flow and 
disturbance determines the amplitude. As in § 2, the disturbance energy- 
balance relationship is required for this purpose. 

If the equations of the mean motion are subtracted from the equations 
(3.1), equations are obtained for the disturbance velocities u’, 7’, w’ (as 
defined by (3.3)). From these it is possible to obtain the following 
disturbance energy equation, which is an exact consequence of the equations 
for an axisymmetric disturbance : 


¢ ‘2 12 2 ei — (ot ra) 
ae To(u “1iq*tiwg 2\r drdz = || | (—pu v It = = \rdrds = 


—p | | (£2 + y'2+4+C)rdrdz, (3.10) 





: 
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the vorticity components of the disturbance being given by 
(rv’). (3.11) 


The three integral terms of (3.10) can be interpreted in a similar way to those 
of (2.5). ‘The net rate of increase of disturbance energy is equal to the 
difference between the integral of the product of the Reynolds stress and 
the flow shear (00/er — d/r), which represents the rate of transfer of kinetic 
energy from the mean flow to the disturbance, and the rate of viscous 
dissipation of kinetic energy. 

The approximate non-linear problem of stability is that of solving the 
equation of mean motion, (3.5), and the equation of energy balance for the 
disturbance, (3.10), where, in the latter, the disturbance is assumed to be 
specified in shape but not inamplitude. The pair of equations (3.5)and (3.10) 
can then be used to determine the amplitude, as in $2. If the amplitude 
of the disturbance in the equilibrium flow is required, @ is given in terms 
of u’v’ by (3.7); then (3.10) involves both second and fourth powers of the 
amplitude, and the latter can therefore be determined. Moreover, an 
assumption similar to that of §2 concerning the growth of the disturbance 
can be made. ‘This assumption is that 0d/dt can be ignored at all times 
on the grounds that it is certainly negligible at both small and large times. 

By the procedure outlined above a differential equation of the form (2.15) 
is obtained, and it has the solution (2.16). In order to calculate the 
coefficients in the equation, it is necessary to specify the shape of the flow 


ou’ , Ou dw’ 1 oa 
eae 
1 


a ~ 
Oz Oz or 


field so that several integrals can be evaluated. In the detailed calculations 
which follow, attention is restricted to the case 2, = 0 so that comparison 
can be made with experiments of ‘Taylor (1936). Secondly, only the very 
simplest case is considered, that in which the annular gap is small compared 
with the radii of the cylinders. ‘he linear stability problem (Chandrasekhar 
1953) is then specified by 

(D2 — )2)(D2 — 2 — o)'v, +2 Tv, = 0, 


(3.12) 
v, = D*v, = D(D*-H-a)o, =0 at C= +}, J 
where 
@ = f4—T; fe = Ty +%s r=r,+Cd, | 
A = ad, a = kd?/v, D= didi, \ (3.13) 
T = OQ; a" y, R= Q), r,d Vv. 


‘The disturbance is proportional to exp(kt). ‘Vhe symbol T denotes the 
‘Taylor number and R denotes the Reynolds number. 
‘lo the same approximation (d <r,) the mean flow is given by 


6 = $7, Q,(1— 22) — (a3 08 d?/v?)O(S), (3.14) 
“f “1/2 “1/2 
Q(¢) = -2 | Z db+2t | Z di+ | Z di, (3.15) 
J —1/2 - -1/2 J -1/2 
Z=-SP, P=(D?-#)S, a, = —ar,Q,P, v0, = (ar? OQ? d/v)S, 
(3.16) 
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a being the amplitude of the velocity u,. The flow shear is given to the 
same approximation by 


06 COD 7,2, 2a?r, Q3d ( ’ : 





( 


or r d y 


Z dt a)). (3.17) 
-1/2 
It has been assumed that the Reynolds-stress term in (3.14) is much larger 
than the terms of order d/r,. 

If the approximation d <r, is adopted and a few transformations are 
applied, the largest terms of (3.10) yield 


da* : y,a we 
n= Ke ys Rea — Se, (3.18) 
where 
7 = 2vt/d?, (3.19) 
"1= | S3dl, y= Z dé, 
J -1/2 —1/2 } 
+ (3.20) 


2 | na 7d) +2 [" ; Pa, y= ; [(D? —d?)P}? dé. | 

This equation is of the form (2.15) and has the property that the amplitude 
tends to a limiting value with the passage of time. If the Reynolds-stress 
term (—y, Ra’) is neglected and the disturbance is neutral, equation (3.18) 
is simply Chandrasekhar’s (1953) variational condition. ‘The critical 
‘Taylor number is given by A?7,, = yq/y. and the equilibrium amplitude by 


4 9 ik 
2 - 2 c 29 
a; 7, (1 7) (3.21) 


Thus a, is proportional to the square-root of the difference between the 
actual ‘Taylor number and the critical ‘Taylor number. 

Equation (3.12) can be used to specity the function S(¢) and, since we 
are primarily concerned with the equilibrium flow, S(¢) will be assumed 
to be the function appropriate to 0 = 0,A =A, = 7, T= T, 1708. ‘This 
is the case of a disturbance at the critical ‘l'aylor number, tor which we 
have, from the variational condition given by Chandrasekhar (1953), 

2 ne it 
a a a a eee ak 


cosh 7f 





~ (3.22) 
cosh bz 


to a good approximation, where A 0-02686. ‘Then for a range of 


‘Vaylor numbers above the critical value, and for the particular wave-number 
\ = 7, the amplitude can be evaluated according to (3.20), (3.21) and (3.22). 


5:425 x 4 708 
g =x “(1 : 7): (3.23) 


This relation, together with (3.16) and (3.22), then gives the velocity 


‘Thus 





distribution of the disturbance. 
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Once the amplitude of the velocity distribution is known, the mean 
skin friction can be calculated from (3.17) and the mean torque on a cylinder 
of given length then follows. ‘Taylor (1936) carried out experiments with 
cylinders 84:4cm long, the outer cylinder having radius r, = 4-05 cm 
while the radius of the inner cylinder varied from case to case. Here we 
shall compare the measured and calculated torque for the case r, = 3-94cm. 
The number of vortices spaced along the axis is 700 approximately and, 
since alternate vortices rotate in opposite directions, the torque measured 
is the mean torque, where the term ‘mean’ is defined in the sense of this 
paper. In figure 5,G denotes the torque measured in units gmcm?sec~ 
and N the angular speed measured in units rev.sec!. It can be seen 
that the agreement is good for T'aylor numbers (T being proportional 
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Figure 5. Comparison of theoretical and experimental torque for flow between 
rotating cylinders. Inner radius 3-94 cm, outer radius 4:05 cm. 


to N?) up to about ten times the critical value (log,, N/v = 1-955), since 
even at log,).V/v = 2:5 the divergence between theory and experiment 
is only about 7",, of the torque. ‘he gradual divergence is probably due 
to upper harmonics of the basic mode of disturbance. ‘The mean velocity, 
calculated from (3.14) is shown in figure 6 for a ‘Taylor number of about 
2:1 x 104. The reversal of sign of the velocity gradient in the centre may 
not be realistic, since this ‘Taylor number is near the limit of validity of 
the theory. The increased velocity gradient at the cylinders is a consequence 
of the fact that the vortex motion requires more power than does laminar 
Couette flow at the same rotational speed. 

From the nature of assumptions made in the present theory, particularly 
the ‘shape’ assumption, it may be expected that the method used above 


B2 
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will be valid tor some range of ‘Taylor number above the critical value, 
and figure 5 shows that this is indeed the case. For a non-linear instability 
problem which has a certain similarity to the present one, that of thermal 
convection in a horizontal layer of fluid, Malkus & Veronis (1958) have 
devised a rigorous expansion procedure which gives the limiting steady 
amplitude of convection for a certain range of Rayleigh number above the 
critical value. ‘Their method gives results very similar to those obtained 
by an integral procedure analogous to the present one. Presumably such 
an expansion would be valid for the cylinder problem also. 
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With regard to the way in which an equilibrium state may be set up 
at a Taylor number above the critical value, it has been suggested here 
(though not proved rigorously) that an infinitesimal disturbance which 
increases exponentially with time will necessarily tend to an equilibrium 
state. On the other hand, the experiments generally require a different 
interpretation concerning growth. ‘Taylor (1923, p.331) states that 
‘The speed of the motor [driving the cylinders] was ... gradually increased 
till instability occurred’? and (p. 342) that “‘ A moderate increase in the 
speed of the apparatus merely increased the vigour of the circulation in the 
vortices...’’. It seems, therefore, that, as the Taylor number ts raised, 
instability sets in at the critical 'Vaylor number with a definite wavelength. 
Further increase of the ‘laylor number entails the flow retaining the same 
wavelength, but with the amplitude of the disturbance motion given by 
(3.23). Consequently, it seems that the non-linear stability problem of 
growth at a given Reynolds or Taylor number, as treated herein, does not 
correspond directly to the circumstances of Taylor’s experiments whereas 
the treatment of the equilibrium flow presumably does. 
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SUMMARY 

The paper deals with the calculation of steady heat transfer 
from a surface of arbitrary temperature distribution to a laminar 
semi-infinite stream of arbitrary velocity distribution. Lighthill’s 
method is improved by a correction which accounts for the 
departures from linearity of the velocity profile within the thermal 
boundary layer, and which comprehends the influences of Prandtl 
number, pressure gradient, body forces, and non-coincident 
start of velocity and thermal layers. Methods are given for 
evaluating the total heat flux directly, and for integrating the 
differential equations for the growth of the boundary layer 
thickness by means of quadratures. 


1. INTRODUCTION 

Approximate procedures for calculating the heat transfer from a body 
to a laminar stream flowing steadily around it fall into two classes. ‘he 
first class contains methods which implicitly assume a fixed relation between 
the thermal and velocity boundary thicknesses. Such is the method of 
Eckert (1942), which has recently been simplitied and extended (Smith & 
Spalding 1958; Spalding & Smith 1958). ‘hese class | procedures are 
valid only for uniform wall temperature. 

Procedures of the second class permit the boundary layer thickness 
ratio to vary; an ordinary ditterential equation is set up for each thickness. 
Such are the methods of Squire (1942), Lighthill (1950), and Schuh (1953). 
Class Il procedures may be used for arbitrary wall-temperature variation. 

The method of Lighthill is asymptotically exact when the thermal 
boundary layer is much thinner than that of the velocity, as occurs, for 
example, when the front part of the body is at the same temperature as the 
stream. Even when the thicknesses are of the same order, the method is 
approximately correct provided that the velocity profile along a normal to 
the surface is nearly linear; this occurs when there is no longitudinal pressure 
gradient. 

Tribus & Klein (1955) have attempted to improve the Lighthill method 
by introducing a correction for pressure gradient due to ‘Tiftord (1951). 
It will be shown that this procedure may actually impair the correctness 
of the calculation. ‘he present paper describes an alternative correction 
which reduces the error of the Lighthill method to 2-5°%, regardless of 


pressure gradient. 
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In addition to this improvement the present paper describes a rapid 
method of integrating the differential equations representing the growth 
of the thermal and velocity boundary-layer thicknesses, and also describes 
a procedure giving the total heat transfer rate. 

It is convenient to consider the case of uniform wall temperature first. 
Cases of variable wall temperature are then dealt with by superposition. 
Finally the calculation of the shear force distribution is discussed. 


2. HEAT TRANSFER FROM A SURFACE OF UNIFORM TEMPERATURE 

2.1. The problem 

Figure 1 illustrates an aerofoil. Measured above the free stream 
temperature, the wall temperature is zero up to a distance € from the leading 
edge, but at larger values of x it has the value 7). ‘Uhe fluid velocity normal 
to the wall is supposed zero, and the fluid properties are taken as uniform. 
Figure 1 also illustrates the growth of the velocity and thermal boundary 
layers. It is seen that the latter is typically thinner than the former, at least 
where x is not much greater than &. 
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Figure 1. Velocity and thermal boundary layers on an aerofoil. 


It is convenient to discuss the properties of the velocity and temperature 
distributions in terms of boundary-layer thicknesses. ‘The nomenclature 
used conforms with that of Smith & Spalding (1958) and is summarized 
in table 1, in which w is the x-component of velocity, u, is the velocity 
outside the boundary layer, T is the fluid temperature, and y is the normal 
distance from the wall. Other notation will be introduced as required. 
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Table 1. 
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The problem is to determine the local and total heat fluxes from the 


wall. It may be considered solved when A, has been determined as a 
function of x; for the local heat flux g” can then be evaluated from the 


- or 1 
gq” k — k a > | 
t ( ey ) 7] 0 A, ( 


where k& is the thermal conductivity of the fluid, assumed uniform. If 
the total heat transfer per unit width of aerofoil, q’, is required, g” can be 


equation 


integrated over the surface. Alternatively, it may be more convenient to 
calculate A, at the trailing edge (~ = /) and then apply the steady-flow 
energy equation in the form 


q = | 9" dx = (cpu, Ty Ae), <1, (2) 


where c is the fluid specific heat, and p is the fluid density. 


2.2. Dimensional analysis 

As is common in the approximate solution of parabolic partial differential 
equations, a ‘ profile’ method is used; that is, the velocity and temperature 
profiles along normals to the wall are assumed to belong to a restricted 
family, in this case those that are found in boundary layers adjacent to 
isothermal wedges in laminar steady flow. For the moment, the shear 
thickness 5, is supposed known as a function of x, and attention is 
concentrated on the conduction thickness A,. 

We assume that the rate of growth of A, in the x-direction depends 
only on the local values of stream velocity, velocity gradient, kinematic 
viscosity v, thermal diffusivity x, and conduction thickness. Standard 
methods of dimensional analysis can be used to restrict the form of the 
relation between these quantities. In performing the analysis it is helpful 
to ascribe different dimensions to lengths in the x- and y-directions; for 
then, by using the boundary-layer assumption that viscosity and conductivity 
are responsive only to gradients in the y-direction, the following continuation 
equation can be directly derived: 


u, dA* ea, A 2 ae 
_ — — , —- , 0 (3) 
v dx v dx 4, 
where o is the Prandtl number v/«. 
This is a first-order differential equation, generally non-linear. ‘The 
function on the right-hand side has yet to be determined. When this has 
been done, standard methods of numerical analysis will yield the dependent 


variable A, as a function of the independent variable x; for u, and 6, are 
already known as functions of x, while v and o are known constants. 
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2.3. The Lighthill method 
Lighthill (1950) has considered the case in which the thermal boundary 
layer is so much thinner than that of velocity that it may be regarded as 
lying wholly within a region of linear velocity profile. ‘To proceed as far 
as possible in the direction of Lighthill’s result by means of dimensional 
analysis, it should be noted that in this case u, and 6, can enter equation (3) 
only in the form of the gradient w,/d,, and that, further, v must have no 
influence. If it is also noted that the application of dimensional analysis 
to the velocity boundary layer gives the continuation equation 
n= 1(3 B). (+) 
where f is an arbitrary function, then the form of the function in (3) can 
be deduced. ‘This deduction leads to the following compact continuation 


equation : 
1/5,\!2 d is \ "7 
BS oa) —| A} —! = const. ( 
x\Uy/ dx O4 


That (5) is a special form of (3), combined with (4), can be verified by 
differentiation of the term within the square brackets. 

By exact solution of the differential energy equation, Lighthill showed 
that the constant in (5) is 6-41. His result was expressed in the equivalent 


mn 
— 


integral form of (5), which, in the present notation, becomes 


1 ] 13] uy 12 fz Uy, 1/2 13 
—_ = —— —— — — 1x i 6 
x7 (ea) we) L1G) 4] o 


(Note that (6-41)!® = (1/3)!3?%, which is the form appearing in Lighthill’s 
paper. ) 


2.4. The Tribus—Klein—Tifford ‘improvement’ 

The Lighthill solution (5) or (6) is asymptotically exact, within the 
validity of the boundary-layer approximation, if A, <6,. It may give 
appreciable error in other cases. Such cases occur among the isothermal 
wedge solutions, when o > 1 and the pressure gradient is finite. Noting 
this, Tribus & Klein (1955), following ‘Tifford (1951), have proposed a 
correction to Lighthill’s formula which depends on o and the pressure 
gradient parameter alone. ‘Their correction is such as to make the formula 
exact for all the isothermal wedges. 

This procedure, however, is open to the following objections. Suppose 
that the heated section on a wedge does not start at the apex but much 
farther back. ‘Then the thermal boundary layer is comparatively thin, 
the Lighthill expression (6) is accurate, and no correction is needed. 
Application of the ‘T'ribus—Klein—Tifford correction will therefore introduce 
an error of the same magnitude as that which it is intended to eliminate. 

It is clear that if the Lighthill formula is to be improved, the correction 
procedure must take account of the relative thicknesses of the boundary 
layers. Yet how is such a correction to be derived from the isothermal 
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wedge solutions if these have a thickness ratio uniquely determined by 
Prandtl number and pressure gradient parameter ?—only by re-grouping 
the isothermal wedge solutions in the light of intuition. 


2.5. The new correction 
Suppose that the correction required, to account for the influence of 
A,/6,, (67/v) du,/dx, and o, is a function solely of the extent of the thermal 
4¢™4 } 1 ’ 
boundary layer into the region where the velocity profile in the boundary 
layer is curved. ‘lhe latter quantity is measured by the dimensionless 
; | : ; 
group (A,5,/v)du,/dx; for the velocity distribution close to the wall can 


be written 
; Gu) (55 . 
( — ‘+ =| — y* 
>) y=0- 2 OV Jy =0 


u, lu,du, , ‘ 
Cia. woke =e (/) 
Og” 2v dx- 


~ 


and at y = A, the ratio of the quadratic to the linear term ts 


(3 qai)/(* =") = 4 Ayo, du (8) 


2v dx * 54 2v dx” 


If the above supposition is correct, it ought to be possible to write the 








continuation equation (3) with high accuracy in the form 


1/5,\'? d (3 _ (- d, du 
—| Ai 6-41+ F( 4 —}], y 
(2) ral 3) | v ae) ) 


in which the form of the left-hand side and the constant on the right-hand 
side derive from the Lighthill solution, and F is the as yet unknown 
correction term. ‘his function may be expected to vanish with its argument. 

The validity of equation (9) can be tested by reference to the isothermal 
wedge solutions ; for if (A,0,/v) du,/dx is as important as has been supposed, 
these solutions plotted in the form suggested by (9) should form a single 
curve instead of a one-parameter family of curves. 

‘This is put to the test in figure 2 on which the isothermal wedge data 
of Eckert (1942) and of Livingood & Donoughe (1955) are plotted. It will 
be seen that the points indeed lie close to a single curve; the supposition is 
apparently well-founded. 

Some scatter is apparent however. ‘The worst is in the neighbourhood 
of (A,d4/v)du,/dx = 0, the case of zero pressure gradient. Here the 
isothermal wedge point for o = 0-7 lies at an ordinate of 6-9 instead of the 
expected 6°41. Since the heat transfer rate is proportional to the reciprocal 
cube root of the ordinate, the error in heat flux is only 2-5°,,, which may be 
considered acceptable. ‘his is of course the same error as was noted by 
Lighthill for the case of the flat plate. 

Use of equation (9) can therefore be expected to yield errors in heat 


flux of less than 2-5", particularly when the heated section starts well 
back from the leading edge. 





i 
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2.6. Free convection and the rotating disc 

Also plotted on figure 2 are points obtained from solutions valid for 
free convection from a flat plate and for forced convection from a rotating 
disc. ‘The former data have been derived from the papers of Schmidt & 
Beckmann (1930), Schuh (1948), Ostrach (1953), and Stewartson & Jones 
(1957). ‘The latter data have been obtained by interpolation in small-scale 
graphs presented by Millsaps & Pohlhausen (1952); their accuracy is 
therefore poor. For both free convection and the rotating disc, the curvature 
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Figure 2. Deduction of the continuation function from the isothermal wedge 
solutions : conduction thickness Ay. 


parameter in the abscissa of figure 2 is —A,(0?u/dy"),, — 9/(eu/ov), — 0, 
instead of the group (A,6,/v)du,/dx which is equivalent for the wedge 
solutions only. In evaluating the rotating-disc solutions, Mangler’s (1948) 
transformation from rotationally-symmetric to plane coordinates was used. 

The free-convection and rotating-disc points lie close to a prolongation 
of the curve for the forced-convection points. Their nearness is further 
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confirmation of the validity of the curvature parameter. ‘The fact that this 
has high positive values for these cases, even with very high Prandtl number, 
is a reflection of the fact that the boundary layers are highly accelerated, in 
the one case by buoyancy forces, in the other by centrifugal forces. 

Figure 2 summarizes almost all the heat transfer solutions available in 
the literature for constant-property boundary-laver flow past an impermeable 
wall. It is evident that more solutions are needed before the field can be 
regarded as adequately covered. Nevertheless, the coordinates used in 
plotting figure 2, by bringing the points near to a single curve, greatly 
diminish the amount of further work which is needed. They may be expected 
to prove helpful also in presenting data for variable properties with a 
permeable wall. 


2.7. Integration 

Many numerical methods are available for integrating equations such 
as (9), in which it will be recalled that A, is the dependent variable, x is the 
independent one, and u,, du,/dx, and 6, are known functions of x. However, 
since F in (9) is normally small, the following quadrature (essentially a 
Picard approximation) will often suffice. ‘The solution of (5) is written as 


3/2 \ 3/2 L , 12 
| 3:(2) | | (5) | - 6412 | (=) dx 
O14, r 04 z 1 
oT, 12 ( \ 
x | (<) F(a, z) dx, (10) 
Je\O, "yp dx 


the subscripts to the brackets denoting the position at which they arc 
evaluated. 

Written explicitly for the required quantity A,(x), (10) becomes (for 
the case when A, = 0 at x = &, as is usual) 


12 12 “4 12 1a sad 
A, (=) | ot x (5!) dx+% | (5) p(A2s z) ax | . (1 | ) 
Uy ¢ O4 Dr O4 V dx 


In (10) and (11), A, appearing in the argument of F is the first approxi- 
mation for A,, obtained from (11) by omitting the second integral. Values 
of F may be read from the curve drawn in figure 2. If in a particular case 
the second integral appears to be important, a second approximation can 
obviously be made. 


2.8. Example 

A calculation of the distribution of heat flux from an ellipse, of axis 
ratio 1:2 and uniform temperature (the case that was first studied by 
Eckert (1942)) has been made by the above method. ‘The result, in non- 
dimensional form, is represented in figure 3, where it may be compared 
with the calculation by Eckert’s class I procedure. ‘lhe present results 
lie somewhat below those of Eckert, but it should not be supposed that they 
are therefore less accurate; for class I procedures involve graver approxi- 
mations than the present class Il method. No experimental results of 


sufficient accuracy are available for comparison. 
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It should be noted that in the example, in which the heated section 
starts at a stagnation point, the integration has to be begun by a straight- 
forward application of L’Hopital’s rule. 
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Figure 3. Calculated local heat flux distribution on ellipse of axis ratio 1: 2. 
3. DETERMINATION OF TOTAL HEAT TRANSFER RATE BY USE OF THE 
ENTHALPY-FLUX THICKNESS A, 

When the local heat transfer has been determined by the method of 
the last section, the total heat transfer rate can be found by integration. 
If only the total heat transfer is required, a quicker method is to evaluate 
u,A, directly and apply equation (2). A continuation equation for A, 
can be set up by a combined use of dimensional analysis, intuition, and the 


exact similar solutions, as before. ‘The solution for the case treated by 
Lighthill may be written 
1 /5,\"" d — a 
oe — (u, A,)?!? = 0-725. (12) 
a\u,/ dx 


The appropriate ‘curvature parameter’ is (A}/?6%/"/v)du,/dx. The wedge 
solutions for various Prandtl numbers have been plotted in the corresponding 
form in figure 4, which confirms expectations by yielding a grouping of 
the points close to a single curve. Fortuitously, this curve is very nearly 
a straight line, so that an approximate linear expression can be found. ‘The 
continuation equation thus becomes 
ies aii 

(3) “Aye 07-01o ee &. (13) 
dx vy dx 

Comparison of (12) and (13) for the case du,/dx = 0 shows the order 
of magnitude of the maximum error involved in the linear approximation, 








30 D. B. Spalding 


and indeed in the method as a whole: about 2:5°, in the solution for heat 
flow. 

Equation (13) can be integrated as a quadrature with a Picard approxi- 
mation as before. ‘he required form is 


a 12 ()- od : \ 2/3 
uA. = | 7 | (3) dx - =i | (u, As)!%5, as | : (14) 
J O4 cf c . 


wherein A, is the first approximation for A,, obtained from (14) by omitting 
the second integral. ‘The curve on figure 4 is so nearly horizontal that a 
second approximation will rarely be worth while. 
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Figure 4. Deduction of the continuation function from the isothermal wedge 
solutions : enthalpy—flux thickness Ag. 


4. HEAT TRANSFER FROM A SURFACE OF VARYING TEMPERATURE 
If the wall temperature 7, varies with € in a known way, the local heat 
transfer at a point xv is determined by superposition in the manner already 
indicated by Lighthill. A brief summary is given here. ‘The effects on 
the heat transfer at x of each increment of wall temperature d7), at & are 
added by the (Stieltjes) integral 
q’(x) =k (“_! grey (15) 
o Ay(é, x) 
where A,(&,.) for each € and the appropriate x is obtained from (7). ‘To 
evaluate (15), x is regarded as fixed in the integration and d7;,(€) is written 
as (d7T,/d&)dé, except at €=( and at discontinuous jumps of wall 
temperature. The wall temperature gradient d7),/d& is supposed given. 
The total heat flow from the aerofoil is evaluated in a similar way from 
the integral 


/ 


q =¢p \, u, A,(é, 1) dTy(€). (16) 
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5. DETERMINATION OF &,/5, AS A FUNCTION OF x 

Finally the distribution of u,/5,, which has been assumed known, will 
be discussed. Tribus & Klein (1955) advocate solving a continuation 
equation for the momentum thickness numerically by the isocline method. 
Simpler methods exist, however: for example, the approximate quadratures 
introduced by Walz (1941), Thwaites (1949), and others. A more exact 
technique similar to that used above for heat transfer will be mentioned 
here. 

Dimensional analysis and the boundary-layer assumption yield a 
continuation equation for the momentum thickness 6, in which the unknown 
function can be obtained from the wedge solutions. A convenient form 
of the equation is found to be 


Ps 83 du 
5-1782) — (). : tipste 7 
ait gy (Mi 35) = 0-4418 —f y= ’ (17) 
for then the second-term on the right-hand side is very small. ‘Table 2 
contains a few values, computed from Hartree’s solution for the wedge 


boundary layers (Hartree 1937). 



































| 
| 62 du, 
—- — 0:0682 0-0266 0 0-0333 | 0:0611 | 0-0855 
Vv dx 
{83 du, ’ . 
t{— 2) — (0-026 —(0-0058 0 0:0033 0-0019 o 
v ax 
| (54/3,)? x 36-9 20-5 |12-75 9-50 7-70 
| 
| | | 
Table 2. 


Equation (17) can be integrated in the manner used for equation (9) 
and (13), a first approximation for 63 being used in the argument of f,. 
A second approximation will hardly ever be worth while. 

Once 5,, the momentum thickness, has been calculated as a function of x, 
the hypothesis that only members of the wedge family of velocity profiles 
occur is invoked for the calculation of corresponding values of the shear 
thickness 6,. For 6,/5, is known for wedges; it is a unique function of the 
pressure gradient parameter, 63(du,/dx)/v. ‘Table 2 gives some values based 
on the work of Hartree. Consultation of a more extended table of 
this nature gives 5, as a function of x, and so enables the heat transfer 
calculation to proceed. 

6, CONCLUSIONS 

1. The Lighthill method of calculating heat transfer in laminar flow has 
been improved by 

(a) a better allowance than that of Tifford for the non-linearity of the 

velocity profile ; 
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(6) the use of the enthalpy—flux thickness for calculating the heat 

transfer over a finite area of the surface; and 

(c) a quadrature procedure for integrating the continuation equation 

for the boundary layer thicknesses. 

2. ‘Vhe effects of pressure gradient and Prandtl number on the rate of 
growth of the thermal boundary layer can be described by a single ‘ curvature 
parameter’ which expresses the extent of the temperature boundary layer 
into the region where the velocity profile curvature is noticeable. 

3. ‘The conclusion 2 is supported by examination of the wedge solutions 


of Eckert and others. 


The author is grateful to Mr D. J. Woodford for working out the ellipse 
example, and to Mr A. G. Smith for many discussions on heat transfer 
through boundary layers. 
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past quasi-cylindrical bodies, with applications to 
wing—body interference 
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SUMMARY 

An extension of the linearized theory of supersonic flow past 
quasi-cylindrical bodies of almost circular cross-section has been 
found which enables a direct calculation to be made of the overall 
forces on wings mounted on such bodies, subject to certain 
restrictions on the plan-form. ‘The method is applied to two 
examples: (i) the effect of an arbitrary body distortion on static 
stability at supersonic speeds; and (11) the effect of wing—body 
interference on rectangular wings mounted on a cylindrical body. 
‘The drag calculations in the second example are compared with the 
results of the supersonic area rule, which is found to be in error 
for moderate values of the ratio of wing chord to body radius, 
though the discrepancy is not serious from a practical point of view. 


1. INTRODUCTION 

‘The linearized theory of supersonic flow past quasi-cylindrical bodies 
of circular cross-section, initially due to Lighthill and Ward (see, for example, 
Ward 1955), has recently been developed by Randall (1955) to cover the 
case when the cross-section of the body is of arbitrary shape differing only 
slightly from a circle. A similar method has also been used by Nielsen 
(1955) in his work on wing—body interference. ‘The method consists briefly 
of expanding the local streamwise slope of the body, together with a suitable 
basic solution of the linearized equation for the velocity potential, as Fourier 
series in terms of the meridian angle; the boundary condition of zero 
normal velocity, which may be considered to apply at a mean circular 
cylinder, is then used to determine the arbitrary coefficients in the series 
for the velocity potential. ‘he Heaviside operational calculus is found 
to be invaluable in simplifying the analysis. ‘The complete pressure field 
due to a body of this type can thus be found, at least in principle, and indeed 
the calculation of the pressure on the surface of the body is now relatively 
simple. But, as Randall (1955) points out, a very much larger amount of 
numerical work is necessary if the pressure off the body surface is required, 
and although much of the basic computation has now been done by 
Mersman (1954) and Nielsen (1957), there is still quite a formidable task 
remaining in any particular case. 
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Such extensive pressure calculations are necessary if it is desired to find 


the effect of the flow round quasi-cylindrical bodies on wings or fins; and 
in effect they are also necessary in problems involving the interference 
potential between wings and bodies. If detailed wing pressure distributions 
are needed, there seems little that can be done by way of further simplifi- 
cation. However, if only the overall forces or moments on the wings are 
required, it has been found that in certain special cases considerable 
inalytical simplification can be achieved and the problem reduced to one 
which is little more difficult or lengthy than that of determining the pressure 
on the body surface. ‘This has been done by integrating the disturbance 
pressure field due to each Fourier component of the body distortion along 
lines in the plane of the wings at right angles to the axis of the body, from the 
body surface to the limit of the disturbance. ‘The restrictions on the wing 
shape which must apply in order that the method shall be valid therefore 
differ according to the particular force or moment in question. In all cases 
the leading edge must be supersonic and the trailing edge at right angles 
to the free stream; if there are subsonic wing tips these must be outside 
the field of influence of the body distortion. Inthe case of drag, it is further 
necessary that the local wing slope should be constant along lines at right 
angles to the body axis, and this effectively limits the method to rectangular 
wings of constant cross-section. 

In spite of these rather severe limitations, the method has direct 
applications to some problems which are of practical interest. First, it 
is shown how to obtain the changes in lift and pitching or rolling moment 
produced by aa arbitrary small distortion of body shape. Secondly, the 
ore conventional problem of wing—body interference is considered, and 
results are obtained for the lift, pitching moment and drag of rectangular 
wings of arbitrary cross-section mounted at incidence on an infinite circular 
cevlinder paralle! to the free stream direction. ‘They are comp ired with the 
previous work of Nielsen & Pitts (1952) and Nielsen (1955), which is shown 
to contain errors for the larger values of the ratio of the wing chord to the 
body radius. 

The results obtained in this way for the wave-drag of combinations 
with rectangular wings and cylindrical bodies are of particular interest 
because they provide an example in which direct comparison may be made 
with the supersonic area rule (Jones 1953). It appears that for moderate 
values of the chord—radius ratio the interference correction to the total 
wave-drag given by the area rule is actually of the wrong sign, and that 
only for very large values of this ratio does it approach the true linearized 
theory. It must, however, be realized that this interference correction is 
quite small compared with the total wing drag (except possibly for wings 
of very small aspect ratio, to which the present method does not apply); 
it is probable that in most problems involving real aeroplane shapes the area 
rule will be much more successful in estimating the total wave-drag, but 


this point still requires further investigation. 
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2. LINEARIZED THEORY OF SUPERSONIC FLOW PAST QUASI-CYLINDRICAL 
BODIES 

We consider supersonic flow with a free stream velocity U past a body 
whose shape differs only slightly from an infinite circular cylinder of radius R,, 
with its axis in the direction of the free stream. We shall use standard 
Cartesian axes Ox, Oy and Oz, with Ow coincident with the axis of the 
cylinder; Oy will later be taken in the plane of the wings and Oz vertically 
upward. ‘lhe origin O will normally be chosen so that the body distortion 
is zero upstream from the plane Ovs (see figure 1). Dimensional coordinates 
will be denoted by (.X, Y, Z); non-dimensional coordinates (x, vy, z) are 
defined by 

x = X/BR,, y= F/R i= Z) Re. 
where 8 = (M*—1)!? and M is the Mach number of the free stream. 
We shall also use polar coordinates (R, 4, -Y) such that 
Y = Rsin#, Z = Reos8, 
and again detine r = R/R,. 

Under the usual limitations of the linearized theory, the motion has a 
velocity potential UN +¢4(X, Y, 7), where the disturbance potential ¢ 
satisfies the equation 
ad ad Ah ad 1 ad 1 a¢ 
axe ay2 oz aR ROR REP 
Introduction of the dimensionless variables defined above reduces this to 


a4 ag a2g a2d ‘a0 1 ag 





(1) 


(2) 


ax2 oy? az? ore ror. r dee: 
We shall use the Heaviside operational calculus (see, for example, 
Jetfreys & Jeffreys (1956) or van der Pol & Bremmer (1950)) in which p 


stands for °/éx, so that equation (2) may be written 


on _7P lop 1 Mp ° dint 
e po fs ee Beer ae i 
pp = r or r a2 ( ) 
(see Ward (1955) for a rigorous justification of this procedure). 


‘The basic solution of equation (3) that is applicable in the present case 


~(p,r, 4) > K.,(pr)ia,(p)cos nt + b,(p)sin né}, (4) 


where K,, is the Bessel function of the second kind with purely imaginary 
argument (as defined in Watson (1944), p.78), and a,(p), 5,(p) are 
arbitrary functions of p. 
In what follows the operational form of a function f(x) will be denoted by f(p), 
where 
f(p) p | e-P* f(x) dx, 
and the relationship between f(x) and f(p) will also be written (cf. van der Pol & 


Bremmer 1950) 


f(p)c f(x) or f(x) Df(p). 
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The boundary condition of zero normal velocity at the surface of the 


body gives 


ad | 


a5 = Un(X, 8), (5) 
OR R R 7 
where 7 is the slope of a meridian section of the body and is assumed to 
be zero for X -. 0. ‘The non-dimensional form of equation (5) is 
1 a 
ae n(x, @). (6) 


UR, oF |, = 


We may expand 7(x,@) as a Fourier series in 0: 


n(x, 0) = > {A,(x)cos nO + B,,(x)sin nO} (x) (7) 


over the range 0 < 6 < 27, where H(x) is the Heaviside unit function. 


The operational form of this is 
(p,8) = > {A,(p)cosné + B,(p)sin 76}. (3) 


n= 
Combining equations (4), (6) and (8) and equating coefficients of cos n/) 
and sin nd, we find that 
a,(p) = UR, A,(p)/pK,,(p), 
and b,(p) = UR) B,(p)/PK,(p). 
Substituting these values in equation (4), we obtain for the final operational 
form of the disturbance potential 
y Slt ‘A (p)cosntb + B (p)sin nd} (9 
~(p,7, 7) UR, > = [A (p)cos ,(p)sin ng}. (9) 
ny PKi(p) 
The function ¢(x,7,@) which is the interpretation of equation (9) satisfies 
the boundary conditions of the body and at infinity, and is identically zero 
for x <r; it is therefore the correct linearized disturbance potential 
over the whole of space exterior to the body. 
The perturbation pressure coefficient (P— Py)/qy (where go = 5po 
py denoting the density of the free stream) is given to the same order of 


Tt? 


approximation by the linearized expression 
Z od 2 
C sae D-ao 
»~ ~ GI?” BR,uP? 
2 K,,(pr) 
5 eee cosnl +B sin nt}. (10 
BD, Kiggy Pathos +B, (p sine 
J 


W(x, r) > W,(p,r) ce’ YB (pr)/K,(p), (11) 
ind we interpret equation (10) by the product theorem (Jettreys & Jeffreys 


1956, $12.13), we obtain 


C, =5 D> 4 cosné | W (x—€&,r)dA (E—r4+1)4 


. sinné | W_(x—€,r) dB, (E—r+1)>, (12) 


ger—I 
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where the Stieltjes notation has been used. Again, equation (12) gives 
the perturbation pressure everywhere exterior to the body. 

Randall (1955) has evaluated the functions W,(x)-= W,(x,1) tor 
r= 1 and x = 0-10 over a wide range of values of x; the determination 
ot the pressure on the surface of the body is thus a comparatively simple 
matter. In cases where W,,(x,r) are required for r > 1, the computational 
problem is much more formidable, though considerable progress has been 
made by Mersman (1954) and Nielsen (1955, 1957). The functions W,,(p, 7) 
detined above (waich differ only by the factor e””” from those used by 
Randall (1955)) are related to the corresponding functions of Nielsen 
(1955) (which will be distinguished here by the additional suffix N) by 


| W,, (é,r)dé =r !2—W,(x,r) for x > 0, (13) 
“0 


and in particular, when r = 1, 


|” W,, dé = 1- W, (x). (14) 


~ 0 


‘lhe results of Mersman and of Nielsen can be used, in conjunction 
with equations (7) and (12), to obtain the complete pressure field of any 
quasi-cylindrical body. 

Unfortunately, it is found that in most cases of practical importance 
(see, for example, Nielson 1955) the Fourier series (12) does not converge 
with sufficient rapidity for small values of (x—r7 +1), so that the pressure 
near the Mach cone x = r—1 1s difficult to determine accurately. This 
fact, together with the large amount of numerical work which still has to 
be done in any given practical case, means that the method is not to be 
recommended as a way of determining overall wing forces if there is any 
reasonable alternative. 

In the next section we shall show that it is possible to obtain comparatively 
simple expressions, which do not in fact require any knowiedge of the 
functions W,,(x,r) except at r = 1, for the integrals of the pressure coefficient 
(as given by equation (10) or (12)) from r = 1 to « along strips perpen- 
dicular to the x-axis. In the following sections some applications to 
practical problems will be given. 


3. DETERMINATION OF CERTAIN PRESSURE INTEGRALS 


3.1. Definition of the functions I,,, ,, 

In problems requiring a knowledge of wing lift and rolling moment, 
it is necessary to determine both the total normal force and the moment 
about Ox produced by the disturbance pressure field acting on an 
infinitesimal strip, parallel to Oy in the plane z = 0, extending from the 
surface of the cylinder R= R, to the Mach lines R = R,+X/f which 
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represent the limit of the disturbance, and of width 6X (see figure 1). 
‘These are given respectively by 





RK \ 


‘ l 1 
éL = OX | (P—P,)dR = —Ry dX | (P —P)\dr, 


and 
“Kt \ -l+er 
6M, OX | R(P— Py) dR = —R,dX | r(P— P,) dr. 
J | 4 1 
since P—P, = 0 for» 1+.x, the upper limit in these expressions may 


be extended to infinity. 


Y 


. 








Figure 1. ‘Typical wing-body combination 


Using equation (10) for (P— P,), we obtain 
| ) 


Rv 'q7'@L/oX D (2/8) > A,(p) | — K,(pr) dr/K,(p), (15) 
r=l 

and Ry “quéM_,/0X D(2/B) & A,(p) | rK,(pr)dr/K,(p), (16) 
where qo | po A | 


hese may be written in the form 


Ro'g7'eL/eX D —(2/B) > pA,(pyl,,,(p) (17) 
and Ry “qi ‘0M _,/eX D3 —(2/B) y pA, (p)L,,.(p), (18) 
where I,,(p) y r'K (pr) dr { pK, (p)}- (19) 


r=] 
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lhe remainder of this section will be devoted to the interpretation of 
the operational functions I, ,(p) 


3.2. Interpretation of I,,, ,(p) 
Lhe transformation pr = € in equation (19) gives 
I... (P) | §"K,(€)db/{p"?K,,(p)}. (20) 
“p 

he upper limit in equation (20) should strictly be poo; but since when 
interpreting I, ,(p) by means of Bromwich’s integral (Jeffreys & Jeffreys 
1956, $12.09) the real part of p is to be taken as positive, and since also 
(A,,€) ~ (7/2€)'"e-* when € > %, the integrals converge at the upper 


limit and are given correctly as written. 
he tunctions | €”A,,(€)d& may be expressed in terms of the associated 


4p 


Bessel functions known as Lommel functions (see Watson 1944, $10.7 


et seq.). We have (Watson 1944, § 10.74) 
| 3"C_(2)dz = (m-+n— 1)2C,,(2)5,,_, »-3(%) — 2C,,_(2)S,, (4), (21) 


where C,,(z) is any cylinder function (arbitrary solution of Bessel’s equation) 
and S$), ,(z) is the Lommel function of the second kind. In particular, 


taking C,(z) to be the first Hankel function H‘’(z), and writing zs = Zé, 
We get 
pet! ) 2 ea (m+n—1EH' (GES, ,,, \0&)—t€H' (2E)S,,,,(1€). 
(22) 
Now A (€) Tl LH (7), (Watson 1944, § 3.7), so that 
| 2K (€) dé [((m+n—1)EK,(E)T,, 1 ,-1(€) + €K,_,(E)T,,,,(E)], (23) 
where Taslth =F?" Ss, (24) 
(he Lommel function S,,,(%) ts a particular integral of the differential 
equation d*y dy 
st tg + (2? —n*)\y = gett, (25) 
dz? dz ; 


When m+n is odd, S,,,,,(%) is defined by the finite series 


f. iP ee lees 1 4 
In i gmk) |] —_— + ill nl 3... |. 26 
When m+n is even, the definition of S,,,(z) is more complicated, but it 
still has the asymptotic expansion (26) when =. 1s large. 


‘ERus: & (3) is a particulat integral of the ditferential equation 


mit 


and 1s given by 
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when m + nis odd, while if m+n is even it has the asymptotic expansion (28) 
when |2/ is large. 
If in equation (23) we make use of the relation 
r _ le an 
— K,_\(€) = K,(€)+n&“K,,(), 
we get 


[enk,(€)dé EK,(E)(mE" T(E) — (m +n — 1)T,,, 4, n-a(€)} 4 
+EK (ET n(€), (29) 


and in particular 
E"K (€) dé = pK, (p)inpT,,,,,.(p) —(m+n— VT, 1,-1(p)}4 
Jp 
+ PK,(P)Tn.n(P), 


since the indefinite integral given by (29) clearly vanishes when & -> 2%. 
Substituting in equation (20), and remembering that W,,(p)= — K,,(p)/K.(p) 
(equation (11) with r = 1), we obtain 


Dan(P) =P" Tn. n(P) —P Wal P) x 
x [np T,,, ,,(p) —(m+n—1)p-T,, 1 »-a(p)]- (30) 
It is convenient to write 
POT (P) = Tnn(P)s (31) 
so that 
Ll P) = Tr( Pp) -—ptW,(p)int,,.,,(p) — (m+ n—1)%, 4p p)}- (32) 


In order to interpret the operational functions T,,,,,(p) and t,,,,,(p) 
we first note that a series expression can be obtained at once from equation 


(28). ‘Thus 
ak (m—1)?—n? 
¥ p? 
{(m— 1)? — n?}{(m — 3)? — n?} 
a ; le 


Tn a(%) — TanP) 





so that 


Tnn(X) = x7/2!+ {(m— 1)? — n?hax4/4! + {(m— 1)? — n?}{(m — 3)? — n?}x8/6! + 


(33) 
If m+n is odd, this series terminates and completely defines 7,, (x). If 
m+n is even, the series (33) still defines 7,,,(%) for sufficiently small x 
(in fact for x < 1); but an alternative approach is preferable, which yields 
remarkably simple results in the cases of interest at present, namely m = — 1, 
0 and 1. 
Consider the differential equation 

( 1) a“ Y _ 2 A (34 

v—l)—54 tas -Wy=- TH, 3 
dx? dx . k! 


where if & is a negative integer the right-hand side is to be replaced by zero. 
We take y to be that solution of equation (34) for x > 0 which satisfies the 
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initial conditions y(0) = A, v(0) = B. Then the operational form of this 
equation may be shown, with the aid of the relation (cf. van der Pol & 
Bremmer 1950, Chs. LV & X) 
x"y(x) D p(—d/dp)"y(p)/p, (35) 
to be 
, ay 


d g 9 
pipe +P gy ~ (P+ my 


ll 


-p k —pB —p*A ifk > | a 
> (36) 


—pB-p*A ifk < 0. 


Comparing equation (36) with the equation (27) satisfied by T,,,,(p), 
we see that they are identical provided that k = —(m+1) and that 
(a) if m= —1, we take A = B=0; (5) if m=0, we take A = 0, B=1; 
(c) if m=1, we take d=1, B=0. Now the further transformation 
v = sing reduces equation (34) to the form 

dy ,  sin*¢ (37) 
— +-nVvr-=—_-=, dI/ 
dd” : k! 
and the solution of equation (37) satisfying the required boundary conditions 
at x = @ = (0) can easily be found. 
In this way we obtain the results 
T_,,(x) = n-*(1—cos nd), | 
7, ,(x) = n'sinn@, > Oc ec! (38) 


T,,(x) = cosnd¢, | 


where tor » = 0 the limiting values as n > are to be taken. 
‘The corresponding results for the functions r,, (x) (equation (31)) 


are 








T_1,(x) = m-*(1—cosnd), (n #0) | ” 
d 
with T_19(%) = 1¢; } (9) 
1 1 1 wg 
T(#) = > oor as Inn 7 + 1)¢ | -1)b (n 1) 
| = , 
with 7),(x) = cosP+sind—1 = (1—x?)!?+xsin'x—-1, 
and = 7) ,(x) = 32°; J 
(40) 
and 7, ,(x) = a cos(n — 2)6— ~ Cos nd — 
vn ee n?>—4 4(n—1)(n—2) 7 2(n? — 1) 
1 
= oo > £08( + 2), (n > 2) 
4(n+1)(n+2)~° 41 
f ) 
with 1, 9(x) = 42°, 
7,(x) = —4+4+ Srsing+ $cos d— 3, c0s 34, 
and = 7, o(%) = 4x? 32+. j 


It may be verified directly that equations (39) to (41) are equivalent 
to the series expressions (33) for 7, ,(%). 
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It remains to consider what happens to 7,, (x) for x 1 in cases where 
r,,,(¥) is not simply a polynomial (i.e. when m+n 1s even) In these cases 
the expressions (39) to (41) appear to break down, but it can be shown 


that they can still be used provided the real part is taken, so that ¢ is to be 


replaced by $7. ‘Vhus when x 1, we have simply 
7, (x)=n?, (n odd) 
(+2) 
and Tighe) = wo; 
To. n(X) = (n*-—1)", (n even + 0) 
(43) 
and Tog *) = bax—1; 
T1,(x) = (n?-—4)"', (nu odd 1) 
(44) 


and T14(%) = f2x 


The range of values of m in equations (39) to (44) can be extended if 
desired by means of the recurrence relation 


a TR 1+ ((m + 1)*—n*},,,,(p), (45) 


whence 


Tin twat) Py? + {(m+1)?—n?! 1 (x—€&)r, (€) dé. (46) 
-~ 
This may be derived either directly from equation (33) or trom the 
recurrence relations satistied by the Lommel functions. 


Returning to equation (32) and interpreting, we get 


Ft) = Fgh) | W(x —&){nz,,,(€) —(m+n—1)ty 1 ,-1(€)} dE. (47) 


Equation (47), together with the expressions (33), (39) to (44) and (46) 
for 7,,,,(¥), enables J, (a) to be calculated without reference to the values 
of W)(x,r) except for 7 l. ‘This has been done tor m= 0 and 1, 
n = 0 to 6 for a range of values of x from 0 to 5; the results are given in 
table 1, and are shown graphically in figure 2. 

In spite of the comparative simplicity of equation (47), it does involve 
certain computational difficulties which are not immediately obvious. 
In the first place, if m+n is odd, 7,,,(v) is a polynomial of degree m+n-+ 1, 
so that except for the smaller values of m+n, /,, ,(x) is given by equation (47) 
as the difference of two tunctions which become large rapidly as x increases, 
while /,,,, itself remains small. ‘This fact, coupled with the oscillatory 
nature of W(x) for the larger values of », makes it difficult to obtain 
satisfactory accuracy in such cases; fortunately, asymptotic formulae ar¢ 
available which are adequate tor most practical purposes (see below). 

‘The second difficulty, which is less serious, arises from the fact that 


when m +n 1s even, 7,,,(¥) has a discontinuity in its first or second derivative 


at x = 1; in tact 7 ,,,(”) becomes infinite at this point. ‘Thus care has to 
be taken in evaluating the integral in equation (47), and formulae of the typx 
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Table 1. 


“0053 


‘The functions J ”" 


i a oe a he ee oe Bs 
TI ae 
0 0 0 O 0 0 
O-O185 0-0177 | 0-0166 | 
0-0675 | 0-0585 | | 0-0457 | 
01055 0-091 0-068 | | 
0-1361 | 0-0999 | | 0-0574 | 
0-2143 0-1257 | | 0-0469 | 
0-3545 | 0-2933 (0-212 | 0-1307 0-068 | 0-0295 | 
0-3664 0-1188 | 0-0219 
0-4293 | 0-0996 | | 0-0261 | 
| 0-6686 | | 0-234 | 0-037 
| 0-4792 | | 0-0822 | | 0-0342 | 
0-5157 | | 0-0721 | 0-0382 | 
0-9958 05396 | 0-194 | 0-0703 | 0-051 | 0-0363 | 
00-5523 | | 0-0743 0-0328 | 
10-5660 | | 0-0805 | | 0-0312 | 
11-3053 0-162 | | 
(05530 | | 0-0856 10-0321 | 
05456 | 10-0881 | | 0-0336 
1-5814 | 0-5358 | 0-161 | 0-0879 | | 0-0344 | 
0-5252 | 0-0860 | 0-0339 | 
05151 | | 0-0838 | 00332 | 
1-8185 | 0-170 | 
0-5064 | 0-0821 | 0-0329 | 
00-4994 10-0815 0-0330 
2-0174 | 0-4944 10-0818 | 0-0333 | 
0-4912 | | 
0)-4897 
0-4895 
00-4904 | | 
0-4917 | 
(-5000 | 0-167 | 0-0833 | 0-050 | 0-0333 | 
(a) 

I, I, » iF Lia | Ke | 
0 Oo | 90 0 ot 6 
01229 | 0-1166 | 0-1067 | 0-094 | 0-080 | 0-065 
0-4709 | 0-3926 | 0-2868 | 0-181 | 0-096 | 0-043 
1-002 | 0-1970 | 0-3718 0-150 | 0-056 | 0-037 

1:6620 | 0-9411 | 0-3596 0-119 | 0-067 | 0-05 
2-4129 | 1-1009 | 0-3303 | 0-127 | 0-070 | 
3-2184 | 1-1932 | 0-3201 | 0-130 | 0-064 | 
4-(527 | 1-2504 | 0-3273 | 0-067 | 
4-8984 | 1-2963 | 0-3340 0-067 
0-3333 | 0-125 | 0-067 | 0-0417 
(h) 


(x); (a) Lyx), (6) L),,(%) 


“ 


] 
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given in Jeffreys & Jettreys (1956, §9.092) must be used near x = 1 in place 
of conventional methods of integration. In spite of these singularities in 
the derivatives of 7,,,,(x), it can be shown that in fact J (x) and all its 





derivatives are continuous functions of x for x 2 0. 
0 eae ae Bi ace — 
5 aaa ] =] ij q ] } 
| | | | | | | 
| | 


45} oe | 4-5 ———__} —__._} | 














(a) (h) 
Figure 2 (a). The functions J), (x). Figure 2(4). The functions J, ,,(x) 


w 


.3. Asymptotic expansions for I,,, (x) 
When x is large, asymptotic expressions for the functions /,, ,(v) can 
be obtained from the definition 


Tes 


I, n(X) D1,,,(p) = —pov-? | em K, (€) dé/K,(p) 
by expanding as a series in p. 
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Now ia n> 0, 
Kale) = M0 DEE *— 21) Wo 2)1(48) 2144 
rannnnrys A +(— 1nd (5g) 

x flog }é+y—4(1+}+...+0-)}+O(E"*loge), (48) 
where y is Euler’s constant. We suppose first that » >m+1. It is 
convenient to denote that part of the series expansion of &”"K,(€) which 
involves only negative powers of € by P(€"K,,). Then €"K,,(€)— P{é"K,, (€)} 
is regular at € = 0, and we can write 


| erK,(@) = | [e"K,(@)—Ple"K Oj] de + | Ple"KO} de - 
p 


Jp 
— |" [enK,,(€)— Per K, (Oj ae. 
“0 
The first term of this series is a constant; and it is easy to see that the sum 
of the last two terms may be written simply as the indefinite integral 
ie 
| €"K,,(€)dé, obtained by integrating the series for €"K,,(€) formally 
term by term. 
‘The result differs slightly accordingly as m+n is even or odd. In the 
former case the important terms of the expansion in ascending powers of p are 
| &€mK,,(£)dé = 2" [(n— 1)! p”—"*1/(m—n +1) + 0dd powers of p] + 
J» 
+ O(p™*" log p). (49a) 


In the latter case a logarithmic term occurs inside the square bracket, due 
to the term in €+ in the expansion of €”K,,(&), and so 


| £mK,(€) dé = 2"[(m—1)!pm—"4Y(m—n +1) + 
- +even powers of p+ O(log p)+O(p"*"*!logp). (49b) 


In both cases 


p*?K(p) = —2" 1[n(n—1)!p-"*! +a power series in p] + O(p” *"*! log p) 
(50) 


Dividing equations (49) by equation (50) and interpreting, we obtain the 


asymptotic expressions 


Fncu(*) ~ 


| 7 ; 
aS | + O(x ~~. if m+n is even; | 
> (51) 
Poy Je ] 1 O xn +1 ; if j is ; | 
ta at | (x )], ifm+nisodd 
When n ~— m+1 the same technique can be used but the pai must 
be obtained individually. ‘The asymptotic expressions for /, (x) and 
/, ,(x) are summarized below; 
Ly o(x%) ~ bax —log 2x—1 + bax + x-*(log 2x — 1/6)— Sax (log 2x— 17/16), } 
i ~ log 2x + x (log 2x — 3/2) + (15/4x")(log 2x — 5/4), f 
I, n(X) — 1; n(n — 1), (n > 1); 


(52) 
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1, (x) = }x*, exactly, } 
I, \(%) ~ 32x — 1 — fax! — 3x2 + Lax 3(log 2x — 19/8), a 
I, (x) ~ | log 2x +! + O(x-), ia 
I, ,(x)~ 1/n(n—2), (n > 2). 

Although these expressions are strictly true asymptotically as x > 2, 


it is found that when ” = | they are unexpectedly inaccurate for moderately 
large values of x, due to certain exponential terms which occur in the 
expressions for W,(x) when n > 0 (cf. Randall 1955). For example, it 
may easilv be shown that 
I,.(P) = 9 : p °W,(p), 
and W(x) contains the term 
2A fa, e* /(1+a7)}, 
where z, ()-6453 +0-50127. We should therefore add to the asymptotic 
expression given in (52) for /) (x) the additional term 
e °6(()-571 cos 0-501x% + 2°174 sin 0-501x), 
and this is found to improve the accuracy considerably. Thus, for x = 3 
the unmodified formula gives /,, = 1°85; the additional term gives 1-53, 
the correct value being 1-58. ‘The corresponding values when x = 4 are 
2:13, 200 and 2-02 respectively. 
In the same way it is found that the correction 
e "1501-797 cos 0-501 4 — 1-756 sin 0-501 x) 


should be added to the expression for /, ,(v) in (53); again this gives an 
accuracy of about 1”,, when x = 4. ‘The corresponding corrections for 


nN 1 are much smaller and may usually be neglected. 


4. 'VHE EFFECT OF BODY DISTORTIONS ON OVERALL WING FORCES AND MOMENTS 

‘The most immediate application of the results of the preceding sect:on 
is to the problem of calculating the effect of small body distortions on the 
overall lift and pitching and rolling moments of the wings of combinations 
for which the body is approximately cylindrical in the neighbourhood of 
its junction with the wings. 

Consider a thin wing mounted centrally at zero incidence in the plane 

() on a quasi-cylindrical body of the type considered previously, such 
that the intersection of the wing leading edges and the body is in the plane 
v = 0; the leading edges must be supersonic, the trailing edges straight 
and parallel to Oy, and the aspect ratio must be sufficiently large that the 


\Iach lines .V 3()°—R,) do not intersect the wing tips (see figure 1). 
We suppose for simplicity that the body is exactly cylindrical and circular 
both for x () and for s — 0, but that the upper half of the body is 
distorted for x () in such a way that the slope of a meridian section at 


the point (.Y,@) is (X,¢). Then only the cosine terms of the Fourier 
expansion (7) need be taken; thus 


n(x,0)= > A,(x)cosnd, (x > 0, 2 > 0), (54) 
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“nt 


Aj(x) = 7! | (x, A)d8, 


where 


“ A, (x) = 271 | 7(x,@)cosné dé, (n> 0). 
‘The disturbance potential for s ~ 0 is then given in operational form 
by equation (9) with the sine terms omitted, since the downwash in the plane 
of the wing is then zero ; for z 0, ¢ is identically zero forward of the 
influence region of the trailing edge since, the leading edges are supersonic. 
‘The pressure on the upper surface of the wing is similarly given by 
equation (10), with 6 = 0 (y > 0) or 6 = a(y < 0). 

Referring now to $3.1, we deduce at once that the /ft on the starboard 


wing is given by 
L | (CL/eX)dX = BRy | (eL/cX) dx, 


~ oO “~ oO 
where C is the root chord and ¢ = C/BR,. ‘Thus the operational form for 
the lift is 


L(c)|(R2q)) D —2 > 1y,,(p)A,(P) (55) 


a] 
from equation (17), in which x has to be replaced by c, so that p now 


corresponds to d/de. 
Interpreting equation (55) by means of the product theorem, we get 


L(c)/(R5 qo) 22 | I,,(c—€)daA,,(€). (56) 
‘The rolling moment .WV, on this wing is similarly given by 


VW (Ri do) 22 I,,(c—&)dA,(€). (57) 


The pitching moment .W/, about Oy is given by 


( 


Vi, | NX(eL eN)dXx PRI | wel eX) dx. (58) 
Now, trom equation (17), 
1 oL 2 e 7 
7 > Io,,(x— €)dA,,(€). 
Ro qo OX Ba}. tole —8)44,(6) 
Substituting in equation (58) and integrating by parts, we get 
V,/(BRido) =2 > E | I, ,(c—€)dA,(E) 


| dx q ae £)44,(6) | 


Reversing the order of integration in the second term, we have 


— 


M,,/(BRiqo) = 2 > | [ely (e—€) — 16, e—€)] dA, (€), (59) 


where I‘~)(€) = |. 1y,(q) ay. 
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he total lift and moment on both wings produced by a body distortion 
of the specified type can now be obtained by carrying out the same process 
for the port wing (taking = 7 in equation (10)), and adding to the 
corresponding result for the starboard wing. ‘Thus 


L\(Req)=-4E[ — Insnle—£)4A,,(€), (60) 
M,,/(Ri Jo) = —4 = . Dion 13(¢— €) dA,,,, i(&), (61) 
n=0 / &=(- 


and 


M, (BR: q) = 4 s | [clo on(C ce £) rar 4 Mc i £)] dA,,(€). (62) 


/ €=0 


5. ‘(HE EFFECT OF WING-BODY INTERFERENCE ON OVERALL WING FORCES 
We shall now consider the application of the theory of §3 to certain 
problems in wing-body interference. In the present paper we shall treat 
only the simplest case, that of a rectangular wing mounted symmetrically 
on an infinite circular cylinder whose axis is in the direction of the free 
stream. ‘lhe basic theory is due to Nielsen (1955), and only a brief summary 
will be given here. 


5.1. Nrelsen’s theory of wing-body interference 
The total disturbance potential ¢ of the combination may be expressed 
in the form 


p oe Pip \ br (63) 


where 6, is the potential of the wing alone (which is usually taken to be 
continued in a suitable manner through the body), :nd ¢; is the additional 
potential necessary to satisfy the boundary conditicus on the body; in 
the present case there is of course no ‘body alone’ potential. 

Provided that the wing leading edge is supersonic, then ahead of the 
influence region of the trailing edge the flows over the upper and lower 
surfaces of the wing are independent, and the interference potential may 
be expanded as a Fourier series similar to equation (4): 


p(p,7,0) = > a,(p)K,(pr)cos nd. (64) 


Here the coefficients a,(p) may differ according as = is greater or less than 
zero, and the coefficients of odd order will vanish when the wing is 
symmetrical about the plane y = 0, as is normally the case. ‘The potential 4, 
defined by equation (63), then automatically satisfies the boundary condition 
in the plane of the wing since 6¢,/¢@ = 0 there. In order that the boundary 
condition of zero normal velocity at the body surface may be satisfied, it 
is necessary that 


Og, /Or = — Opy/or = — Rou,y, (65) 
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when r= 1. Here w,, is the radial component of velocity on the surface 
of the body due to the wing alone. Let u,,, be expanded as the Fourier 


series 
ig = =U y A,,,(x)cos 2n8, (66) 
where a 
Ay = —7 1 " idea U 1 dé, 
and A. = —2n [" u,y U-lcos 2nd, (n> 0). 


-~ UO 


Then, if A,,,(p) are the operational forms of A,,,(x), the boundary condition 
(65) gives 


>,(P) = To l 'A,,(p) (PKz,(P)}, (67) 


so that the interference pressure field in the plane of the wing is given by 


Cyr = (Pr Po)/do = —28* & Aa(P)Ken(Pr)/Kon(p)- (68) 
n=0 

It is clear that the ratio —u,,,/U corresponds to the body slope 7 
considered in §2, $3 & $4, so that the effect of the interference between the 
flow field of the wing and the body is equivalent (so far as the wing is 
concerned) to a simple distortion of the body of the type considered 
previously. ‘The remainder of the theory is thus exactly as given in §2 
and need not be discussed further here; some results concerning the 
detailed pressure distribution are given by Nielsen (1955). 


5.2. Lift and pitching moment of rectangular wings on a cylindrical body 
Consider a thin rectangular wing, whose aspect ratio is sufficiently 
large that the Mach lines from the leading edge of the wing—body junction 
do not intersect the wing tips (the aspect ratio A of the exposed wing must 
exceed 2/8). ‘The wing is mounted at incidence « on a cylinder of radius ry 
which is at zero incidence, so that the wing leading edge coincides with the 
y-axis. ‘The flow field of the wing alone is then two-dimensional in the 
appropriate region of the body, and consists simply of a constant downwash 
U over the region zx < z. We have therefore u,. = —Uzsin@ for 
z x, while wu, = 0 elsewhere; and the Fourier coefficients A,, (x) 


are given, for 0 <@ a, by 


A, = 2471 | sin@ dé 


- A,, = 4an7 | sin#cos2n6 dé, (n > 1) 
where ¢ = sin-!w if x <1,and¢d=}z7ifx>1. Thus (see Nielsen 1955) 
Ay, (%) = af,(%), 
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where 
fox) = (2/m){1-—(1—x?)"*} if x <1, 
2) ft “> 1 
and, for n is \ (69) 
. 2fcos(2n—1)6 cos(2n+1)¢ 2 fx <l 
fau(%) 7 2n—1 ~  On+] “13 * : 
4.n(4n?—1) — on 2 > ih. 
§.2.1. Lift 


The total lift increment due to interference, which is twice that on the 
upper surface of the wing, is given by equation (60) as 
L,(c) (Rio qo) 2855 Io.2,(¢— €) dAz,(€). (70) 


a 
uv & 0 


Since A,,(x) are all continuous functions of x such that A,,(0) = 0, 


equation (70) gives 


L,/(Rigo%) = —8 > | Ine,(e—€)fo,(€) d€, (71a) 
where 
f(x) = 2m-1x(1 — x?)-12 
re (x) = 40-1x(1—2x?)-!*2 cos2nd (n > 1) | 
when 0 < x 1, and f,,(x) =0 for x >1. Thus, if c >1, 
L,|(R2qo%) = —8 > | Igan(—fan(2) e. (71b) 
We may write this result 
L; (Ri qo) = —8G(c), 
where 
He= > Gols, 
and ipa : (72) 


[t will be noticed that the functions /f;,,(~) become infinite like (1 —x)-!2 
when x = 1; because of this singularity the numerical integrations up to 
€ = 1 were completed by using a formula of the type given by Jeffreys & 


Jeffreys (1956, $ 9.092): 


| (ax? + Bul? + yx3? + 6x°?) dx = h(5-994y, — 8-S36y, + 8-974y, — 


—1-854y,), (73) 
where y,, is the value of the integrand at w = nh. 


The functions Gy, G,, Gy and G, have been computed for a range of « 


from 0 to 5, and are given in table 2. 
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al Go Gs | G, } Gs | 


Pe Pe SE te 
0} 0 o | o | o | 
}0-2} 0 0-0001} 0-0001} 0-0001 


| 0-4 | 0-0006 |} 0-0012} 0-0009} 0-0005 
| 0-6 | 0-0032 0-0052 0-0023 | -0-0003 
‘8 | 0-0102 0-0135 0-0008 0-0049 
0 | 0-0258 0-0234 | —0-0098 | —0-0065 
-2 | 0:0583 0-0214 | —0-0164 0-0027 
‘4 | 0-1088 0-0015 0:0165 0-0010 
6 | 0:1743 0:0294 0-0133 0-0015 
*8 | 0-2528 0-0663 0-0087 0 

‘0 | 0:3424 | —0-1046 0-0052 -0-0020 
-2 | 0-4417 | —0-1407 | —0-0036 | —0-0026 
-4 | 0-5494 0:1723 0-0039 0:0017 
0-6645 0:1979 0-0053 | —0-0008 | 
0-7861 0-2170 0-0068 | —0-0008 | 
‘0 | 0-9135 0-2298 0-0079 | 


CQ 











— 





Fe WWW WWD NN WN ND 








2 | 1-0461 0-2371 0-0082 | | 
4 | 1-1833 | —0-2397 | —0-0081 | 
| 1-3245 | —0-2388 | —0-0076 | | 
8 | 1-4695 | —0-2355 | —0-0071 | 
0 | 1-6179 | —0-2309 | —0-0068 | 

4-2 | 1-7694 | —0-2258 | 

4-4 | 1-9234 0-2208 | 

4-6 | 2-0801 | —0-2163 | 
4-8 | 2-2390 | —0-2127 | | 
| 5-0 | 2-4000 | —0-2100 | | 

| 2122 | —0:-0071 | 0-0012 


| oO | 0: 





Table 2. The functions G,,(c). 


When c is large, an asymptotic expansion for G(c) may easily be obtained 
from the results given in equation (52) for J) ,,(x). We have, when ¢ > 1, 
integrating equation (72) by parts, 


“1 ) 
Gy(c) = (2/7) sLoo()— | Loole—€)(1 - )-* dé >, 


and since Ip (x) ~ $ax—log2x—1-+ }x/x + O(log x/x?), it follows that 


G,(c) ~ c—(2/z)log 2c — (2/7 + 4a) +c + O(c log c). (74) 
Similarly, we find that 
G,,(c) ~ —2/n7(4n? — 1)(2n—1) (75) 
when n > 1; and since it can be shown that 


¥ [n(4n? —1)(2n—1)}1 = 1x? +} —2log2, 


it follows that 
G(c) ~ c—(2/z)log $¢e— (3/7 + $7) +71 + O(c * log c). (76) 


It is found that for c = 5 the asymptotic formula (76) gives a result which 
differs by less than 5°,, from that obtained by direct computation, so that 


D2 
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it provides a satisfactory method of extending the calculations to any desired 
value of c. Greater difficulty is experienced for small values of c, due to 
slow convergence of the series }G,,(c). A first approximation as c 
approaches zero can be found, following Nielsen (1955) by considering 
the body as a vertical boundary on which is a given distribution of sources 
corresponding to the normal velocity produced by the wing. In this way 
we find that 

G(c) = &/9x when c > 0. (77) 


The use of four terms only of the series ¥ G,,,(c) in determining the function 


G(c) appears to be adequate, with a probable error not exceeding about 2°, 
T T ois | aa | 
A 
/ 
/ 
4 1 
/ 
/ 
/ 
/ 
/ 
/ 
/ 
J 
4 
if 
Figure 3. The function G(c). 
for « 2; for smatier values of c the accuracy decreases rapidly, while 


the formula (96) cannot be expected to apply for c > 3. ‘There is thus an 
intermediate range in which G(c) has to be estimated graphically, but this 


is not expected to introduce any serious error. The function G(c) is shown 


in figure 3. 
Nielsen (1955) expresses his results for the effect of interference on the 
total lift of rectangular wings of finite aspect ratio in the form of a coefficient 


k,,, defined as the ratio of the lift on the exposed wings in combination 
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with the body to that on the exposed wings joined together in the absence 
of the body. It can easily be shown that, provided the exposed aspect 
ratio A exceeds 2/8, so that tip effects do not enter into the interference 
calculations, the ratio k,, may be expressed in terms of the function G(c) 
- Ryy = 1-—4G(c)/c?(2BA — 1). (78) 
The quantity +G(c)/3c?, which gives the value of (1—,,;) when BA = 2, 
is shown in figure + and compared with the corresponding result obtained 
by Nielsen (1955). ‘The agreement between the two methods is reasonably 
good for the smaller values of ¢ (less than about 4), but Nielsen has over- 
estimated the magnitude of the interference effect when c is large. ‘This 
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Figure 4. Effect of wing—body interference on the lift of a rectangular wing on a 
cylindrical body. 


is due to the use of an incorrect asymptotic formula, which is equivalent to 

G(c) ~ (2c/7)(1 + log 2) — (2/z)log ¢ + O(1). (79) 
The error appears to be due to the fact that equation (79) was obtained by 
integrating an asymptotic expression for the span loading; and though the 
latter is correct, an examination of its derivation shows that it is applicable 
only when x/r (or more precisely (v—r+1)/r) is large. Thus, even when 
c is large, this expression is only correct in a region near the body, and it 
is not permissible to integrate it in a spanwise direction to obtain the total 
lift. 


5.2.2. Pitching moment and centre of pressure 
The increase M,,, in nose-up pitching moment about the leading edge 
due to interference is found by means of equation (62) to be 


D ec 


My1|(BR 40%) = 8 p | [Toen(e—€)—Lh 2 (¢—€)]fon(é) dé, (80) 
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and this is equivalent to 
M,, (BR q)%) = 8[ceG(c) —G~(c)], (81) 

where — 

G—“o)\-= | (Glee)dx: 

Using equation (72) for the lift decrement, combined with the standard 
linearized theory for wings of finite span, it can easily be shown that the 
position of the centre of pressure is given in terms of the chord by 

Xo} l | - + 12c-*G(c) — 24c 3G aa 


6k, y(28A 1) 


where k,y is given by (78). ‘The asymptotic expression when c is large 1s 


_——- 1 + 24e-(log 4c + 4? — 3) e 
a i —~“e i-CF — 


‘The variation of the position of the centre of pressure with ¢ is shown in 
I 


(82) 





figure 5 for BA = 2; again Nielsen overestimates the forward movement 
for c greater than about 4. 











Figure 5. Effect of wing—body interference on position of centre of pressure for a 
rectangular wing on a cylindrical body. 


5.3. The drag at zero lift of rectangular wings on a cylindrical body 
‘The theory of $5.1 can easily be extended to give the drag of any 
rectangular wing ot symmetrical cross-section at zero incidence. We suppose 
that the shape of the wing section is given by the slope 
dZ/dX = &(X/C) for z > 0. (84) 
The wing alone produces, for z > 0, an upwash field 
u.y = Ud(X/C—BZ/C). 
(No connection is intended with the Dirac delta function; 6(v) may be 
an arbitrary function of x in 0 <x < 1 and may have a finite number of 
discontinuities in that range.) This upwash field has a component normal 
to the body surface : 
uy = U sind d(x/c— sin /c) when-d <d0<¢ 
andz—d <@0<7+4d, > (85) 
- 0 elsewhere, 


where as before 4 = sin’! x if 0 v 1 


and ¢ = }rifx> 1. 
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The Fourier components of —u,,,/U are thus 
2 /¢ . ,./x—sin6\ 
A, = —- | sin06 ~~ dé, 
TJ c 
; (86) 
and 4 /? . aofX*—sind 
‘ A,,=-- | cos 2nd sin 86" dé. 
2 =a : 


If we write sin@ = € and cos 2n6 = C,,,(€)(a polynomial of degree 2n in &), 


4px ae 
A, (x) — es | yrs Col 8(* , *) dé. 


sut (cf. equation (69)) 


then 








fo, (x) =- | sin@cos2né dé, 
47. Culé) 
ss }§ ian’ 
Hence 
’ x—€ 
hie = | fa(O2( - ) ae 
— fe,,(%)d(O) — - | fon(x — €)8'(E/c) d&. (57a) 
[his may be written in Stieltjes form 
A,,,(x) ak fo, (x — €) d&(E/c), (87 b) 
where 4(€) is defined to be zero for € < 0. 
The total interference drag D, is given by 
-C ~ o 
D,=4| &X/C)aX | (P,—P,)dR 
x=0 / R=R, 
2R* py tie d(x/c) dx | BC,,,(x, r) dr. (35) 
JZ 0 r=] 
Now we consider the operational form 
| BCyrdr=2 ¥ pAs,(p)o2,(P) (89) 
2 Ld 5 \ I), o,,(x — €)A5,, (€) dé (90) 
eh OU ee 


interpreting equation (89) by the product theorem and making use of the 
fact that A,,(0) =. Substituting in equation (88) and integrating by 
parts, we get 

9 , dz 
D,/(R5q) = 8 | A(x c) dx = | 


1x 2 fo; 
0 “0 


y 


n(x — €)Ap, (€) dé 


-7 


=-8| d&(x/c)|° > Ino,(x—8)A44,(E)dé, (91) 


/ r=0 “0 n=0 


where again 6(&) is defined to be zero for € > 1. 
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But, from (87), 


[ Tyo, (x — €) Az, (E) dé ; Toon(x — &) dé f, (€ — ©) dd(C/c) 
~{ — d&(f/c) | Iya, x — Ef, (E- dé 


— | Gy,,(x — €) d5(E/c), (92) 


where G,,,(x) is defined above (equation (72)). Substituting in equation (91), 
we have finally 
D,/(Roqo) = 8 | dd(x/c) | G(x — €) dd(E/c) (93 a) 
I 0 - é 0 
It is necessary to use the Stieltjes form (93 a) if 6(€) has any discontinuities 
in the range 0 — & 1, but if 6(€) is continuous in this range, then the 
extended form of this equation is 
D,|(8R5 qo) = 67 | 8'(x/c)dx | 8'(€/c)G(x«— &) dé + 
“0 “0 


te} 6(0) | 8/(x/c)G(x) dx — 6(1) i 6'(1 — x/c)G(x) dx 


~ 0 ~ 0 


d(0)0(1)G(c) Pe (93 b) 


For a section symmetrical about x/c = }, this simplifies to 
D,/(8Ri qo) = 7 | 8'(x/c) dx | 8(E/c)G(x— €) dé + 


0 0 


+ 2¢ 18(0) | 8’(x/c)G(x) dx + 6°(0)G(c). (93 c) 
For a wing of biconvex parabolic section and thickness ratio 7, we have 
8(x/c) = 27(1 —2x/c), (94) 


and equation (93 c) gives 


D, (R> do) 877, 16c-2 | dx | G(E)d&—16c-! | G(x) dx+4G(c)> 


0 ~ oO 


872. 1l6c-? | (c—x)G(x)dx—16c-! | G(x) dx +4Gi(c)> 


0 ( 


-32c-*774 4 | x G(x) dx —c?G(c) >. (95) 


0 


The interference drag ratio k,,,, defined in a similar way to ky by 








P Dy T D, Total wave drag 
_ Dd. Wave drag of exposed wing alone’ 
is thus given by 
6 4 ¢* : : 
Rny 1 — BA 3 } vG(x) dx Gi(c) . (96) 
BAc* | c* J 


tor a wing of biconvex section. 
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The corresponding result for a symmetrical double wedge section (for 
which the Stieltjes form (93a) has to be used) is 


4 
Rpy = 1- BAe (4G(3c) — G(c)}. (97) 


k 


These results are shown in figure 6, where the drag ratio kp, is plotted 
against the equivalent chord-radius ratio c = C/(Sry), for wings of exposed 
aspect ratio 2/8. As ¢ increases from zero, kp, at first rises from 1-0 to 
a maximum of about 1-03, reached when c is between 3 and 4; it then 
decreases steadily, being again equal to 1-0 for c between 6 and 7, and 
reaches a minimum of about 0-97 when c is 15, after which it increases 
slowly towards the asymptotic value 1-0. The difference between the 
results for the two wing sections is very small. 
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Figure 6. Effect of wing—body interference on wave drag of a rectangular wing on a 
cylindrical body. 


The results for a double wedge section may be compared with those 
obtained by Nielsen (1955); there is good qualitative agreement but again 
Nielsen has overestimated the magnitude of the interference effect, for the 
same reason that has been suggested in the case of lift. 

A qualitative explanation of the behaviour of the interference drag 
can most easily be given by considering a wing of double wedge section. 
The forward half of the wing produces a compression wave which is in 
effect reflected from the bedy as an expansion acting over the influence 
region of the junction of the wing leading edge with the body. ‘This 
expansion causes a thrust on the forward half wing and a drag on the rear 
half. Similarly, the rear half wing produces an expansion which is reflected 
as a compression over the smaller influence region of the junction of the 
half-chord line with the body; this produces a further thrust. ‘The balance 
between these components of thrust and drag depends on the value of c; 
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and it may be seen from consideration of the corresponding results for a 
single wedge (which are qualitatively similar to those for the lift ratio Rk, yy 
(figure 4+)) that for small values of c the drag on the rear half outweighs 
the two components of thrust, while as c is increased the balance is gradually 
altered until for large values of c there is a resultant decrease in the total 
drag as compared with that of the wing alone. ‘This balance between 
opposing forces of thrust and drag also explains why the magnitude of the 
effect of wing-body interference on drag is in this case very much less than 
on lift; thus the drag correction does not exceed + soe. while the lift 
correction has a maximum value of about 12°, (both for BA = 2). 

An interesting general result for large values of ¢ may be obtained by 
substituting the asymptotic expression (95) directly into the interference 
drag formulae (93). It may easily be shown that the leading terms of the 
resultant asymptotic expression for D, are 

1 
D,/8Rigy ~ —c | [8(€)]? d€ +7! log c(0?(0) + 67(1)} +O), (98) 
0 
provided that 6(€) is continuous in 0) € iP 
Now the wave drag of the exposed wing panels joined together is given by 
1 
Dy ldo = 48 *AC* | [OCE)) de, 
~ O 
l 
Dyy/(8R5 qo) = BAc® | [8(E)]*? dé. 


so that 


Hence 
: b ioe Oe 


; | 
b oe oO 
ha Oo?g 


+ O(c!) >. (99) 


It is noteworthy that the first term in this expression is independent of the 
section shape. For a biconvex section, equation (99) gives 


Rnyw ~ —2 BAc\1— (6/zc)log c+ O(e~')}, (100) 


which may be verified from equation (96). 

It may also be shown from equation (97) that the first two terms of the 
corresponding asymptotic expression for k,,, for a double wedge section 
are in fact identical with (100), but it appears from equation (99) that this 
is not true for a general section shape. 

Although it has been shown here that the effect on the total wave drag 
at zero lift of the interference between rectangular wings of moderate aspect 
ratio and cylindrical bodies is small from a practical point of view, the 
problem considered here is of some fundamental interest because it provides 
an example in which it is possible to compare the results of true linear theory 


with those of the supersonic area rule. 


6. COMPARISON WITH THE SUPERSONIC AREA RULE 


It has recently been realized that the supersonic area rule due to Jones 
(1953) for the wave drag of wing—body combinations is not strictly correct, 
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even within the limitations of the linearized theory, except for a very 
restricted class of objects, which it is not at present possible to define 
at all clearly. ‘This question has been discussed in detail by Lomax & 
Heaslet (1956), who give several examples in which the area rule can be 
shown to lead to considerable errors; the problem considered in § 5 provides 
a further example. 

The area rule is in effect based on two fundamental assumptions, which 
are to some extent inter-related: (a) that the flow round a symmetrical 
wing-body combination at zero lift may be represented by a distribution 
of simple sources over the plane of the wings and along the axis of the body; 
and (6) that the total (linearized) velocity potential for the combination is 
equal to the sum of the potentials for the exposed wing and the body taken 
separately, each in the absence of the other. 

‘These assumptions are clearly incorrect in the case of wings mounted 
on a cylindrical body. For such combinations assumption (4) implies 
that the velocity potential is simply equal to that of the exposed wings alone ; 
and this will not in general satisfy the condition of zero normal velocity on 
the body, so that an additional axial distribution of poles and multipoles 
has to be added, thus violating assumption (a). 

In order to make a more direct comparison with Nielsen’s method, the 
area-rule potential may be thought of as made up of the potential of the 
complete wing continued through the body, together with the negative 
of the potential due to the portion of the wing blanketed by the body; 
and the latter should be equal to the interference potential 4, of $5.1. 
Although these potentials are to some extent similar, it is obviously 
impossible in general for them to be completely equivalent. ‘Uhat this is 
so in the case of rectangular wings may be verified indirectly by comparison 
of the interference drags obtained by the two methods. 

Consider a rectangular wing of overall span 2s and unit chord, mounted 
on a cylindrical body of radius Ry. The oblique area distribution S(é, 4) of 
the exposed wing, defined as the projection on planes perpendicular to Ox 
of the area of a section by the plane 

X = £+BYcos8, (101) 
is clearly given by 
S(E, 0) = S(E, 9, s) — S(E, @, Ro), (102) 
where S(é,6,s) is the corresponding area distribution for a complete 
rectangular wing of span 2s. 
The derivative of S(E,6,s) with respect to € is given by 
S'(E, 0,5) = 281 sec Of f(E + Bs cos 6) — f(E — Bs cos 4)}, (103) 


where Z = f(X) is the equation of the wing section (cf. Lock 1957), and 
{(X) is defined to be zero outside the wing chord. 
The supersonic area rule gives for the overall wave drag 


Diq, = -—7 l" dé 1 [ dédn log €—7 S"(E, @)S’(, 4), (104) 


0 
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where primes denote differentiation with respect to € or 7; the double 


integral is taken over the range of € (or 7) for which S’(€, 9) does not vanish. 
Substituting for S(é,@) from equation (102), we see that 


D = D(s)+ D(R,)—2DAs, Ro), (105) 
where D(s) is the drag of a rectangular wing of span 2s with the given section, 
and D; is defined by 

qi Ds, R,) 7T : | 4 dé | dédy log e - 7) © te; GO, s)S"(m, 6, Ry) (106) 


—47-°8-* | sec?6dé | | dé dnlog €—n 
’ {f'(E + G;) -f'(E-a) HI (a + o2)—f ( — O>)}; (107) 
where o, = Bscos@ and o, = BR,cos#. Linear transformations of the 
variables € and 7 reduce this to the form 


n/2 ‘1 fl 
qo 'DAs, Ro) = —47°282 | sectadd | | dx dy f’(x)f'(y) » 
“ ) 


~ 0 ~ f 
x log! {(x — y)? — (a, — a9)? }/{(w —y)? —(0,+5)?}; (108) 
and by comparing this with the corresponding expression for D(s) (see 
Lock 1957) it may easily be shown that 
D; = Dii(s + Ro)} — DUL(s— Ro)}. (109) 
Now if s > }8-', it is known that 
D(s)/qy = 2sCp,, 
where C,, is the two-dimensional wave drag coefficient for the given wing 
section; andif the aspect ratio d = 2(s— R,) of the exposed wing is restricted, 
as in $5.3, to be greater than 2/8, then equation (109) gives simply 
Dilgg = 2RoC 5. 


Substituting in equation (124), we get 








D/ qo = 2(s — Ro) Cp. + {g-'D( Ro) — 2Ro Cp,}- (110) 
‘The first term represents the drag of the exposed wing panels joined 
together, so that the interference drag ratio k,,, is obtained as 
ZR, Co 2 Ga, 
knw = 14 a ( 1 l+-s (= “4 (111) 
(Ss R,) ( D BAc C D 
where C, is the drag coefficient of the portion of the wing blanketed by 
the body, whose aspect ratio is 2/(8c) in the notation of $5. (This result 
has also been obtained by Sheppard (1957).) Equation (111) suggests { 
that kp, = 1 for 0 <c <2, and that R,, 1 for c>2. This can 


easily be seen qualitatively by considering the physical assumption on 
which the area rule is based: that the flow field is equivalent to that of the 
exposed wing panels with the body removed. When c is less than 2, there 


is no interaction between the two wings and the drag is unaltered; as soon 
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as c becomes greater than 2, the compression wave from the inner leading 
edge of one wing acts near the trailing edge of the other to give a reduction 
in total drag. 

The value given by equation (111) for kp, according to the area rule 
has been calculated for wings of biconvex and double wedge section, using 
the results of Harmon (1947); we obtain 


2 zs } 4 4\12 
Rpw = 1- oF el ne - “| (6 = )eosh Me-(1 - 3) ] (112) 


for a biconvex section, and 
2 
py = 1- npAc (ce) -f(3e)}, 
TOl1¢ 
where 


f(c) = W(c2 —4)"2 + (2/c)cosh Le — 2 cos“1(2/c) (113) 


for a symmetrical double wedge section. 

These results are included in figure 6 for comparison with those of the 
accurate linearized theory of §5. It is clear that for values of ¢ less than 
about 6, the area rule is completely erroneous in predicting the effect of 
wing—body interference in the present case, and that not until c exceeds 10 
is there any quantitative agreement with linearized theory. ‘The latter 
point was to be expected from previous work on the subject; it has frequently 
been stated, and recently proved rigorously by Fraenkel (1958) that the area 
rule is always true for combinations of the type considered here, provided 
that the ratio corresponding to c of the present paper is large. In fact it 
may easily be verified that the first two terms of the asymptotic expansions 
for kpyy when ¢ is large, obtained from equations (112) and (113), are 
identical with the corresponding terms in equations (96) and (97) (cf. 
equations (100) et seq.). 

It must again be emphasized that the actual numerical difference between 
the results of the area rule and of exact linearized theory is not large for the 
example considered here; the ratio Drag (Linearized Theory)/Drag (Area 
Rule) has a maximum of about 1+0-12/84 when c = 4, for values (> 2) 
of BA which are covered by the present theory. But there is little doubt 
that for smaller aspect ratios the error in estimating the wave drag of the 
wings by the area rule will continue to increase and may become of practical 


importance. 


7. CONCLUSIONS 

An extension has been developed of the linearized theory of Nielsen 
(1955) and Randall (1955) which enables more rapid and accurate calculations 
than hitherto possible to be made of the overall forces and moments acting 
on a restricted class of wings mounted on quasi-cylindrical bodies. ‘Two 
principal examples have been studied in detail. 

The first of these is concerned with the changes produced by small 
distortions of an otherwise cylindrical body on the static stability of thin 
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centrally mounted wings at zero incidence, subject to the condition that 
the trailing edge is straight and unswept and that the leading edges are 
supersonic. 

The second application is to the problem of wing—body interference at 
supersonic speeds, in particular to the example previously considered by 
Nielsen (1955) of a rectangular wing mounted on a long cylindrical body 
of circular cross-section, the latter being at zero incidence; both methods 
are applicable only if the exposed ratio exceeds 2/8. ‘The ettect of 
interference on lift and pitching moment ts found to be in reasonably good 
agreement with Nielsen’s results for small values of the equivalent 
chord-radius ratio « C/8R,, but for values of c greater than about 4 
the theory of the present paper, which predicts smaller interference etfects 
than does Nielsen’s (due to a slight error in the latter), should be the more 
accurate. 

A new formula has been derived for the ettect of interference on wave-drag 
at zero incidence which enables this to be computed rapidly and accurately 
for a rectangular wing of arbitrary cross-section. ‘The results for double 
wedge and biconvex parabolic sections are compared with those of the 
supersonic area rule, and it is found that there is poor agreement until the 
ratio c is greater than about 10; and this is a value that is seldom attained 
in practice. ‘Che numerical discrepancy in the total wave drag of the wings 
is not large for aspect ratios covered by the present method, but may well 
become more serious for values of BA much less than 2; although, since 
only the drag of the wings has been considered here, the area rule may still 
give reasonable accuracy for a typical complete wing—body combination 
when the drag of the fore- and after-bodies is taken into account. 

It is clear that further comparisons of the type given above are desirable 
before a true assessment of the accuracy of the supersonic area rule can be 
made. One possible extension of the methods of the present paper would 
be to study the effects of body modifications on the wave drag of combinations 
with rect ingul: ir Wings, and it may also be possible to design body shapes 
to give an overall reduction in drag; these shapes and the corresponding 
drag values could then be compared with the predictions of the area rule. 


‘he work described above was carried out in the Aerodynamics Division 
of the National Physical Laboratory and is published by permission of the 
Director of the Laboratory. 
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SUMMARY 

The initial value problem presented by mixing and chemical 
reaction in the wake of a flat plate is solved using the boundary- 
layer approximation. When a cool combustible mixture and its 
hot combustion products are separated by a finite, perfectly 
insulating flat plate, the velocity, temperature, and combustible 
concentration are determined in the vicinity of the trailing edge. 

The mixing problem without chemical reaction is solved in 
terms of a ‘universal solution’ for a given initial temperature 
ratio and Prandtl number from which the solution for arbitrary 
temperature ratios can be obtained. 

The mixing problem with chemical reaction is solved in terms 
of a ‘universal solution’ for the first two terms of an assumed 
series solution for the temperature. In this case the ‘ universality ’ 
is with respect to a parameter B characterizing the chemical and 
hydrodynamic initial conditions. 

The axial distance from the trailing edge to the first local 
temperature maximum is given in terms of the initial conditions 
and is shown to be greatly shortened by the presence of the viscous 


wake as compared with non-viscous mixing. 


i. INTRODUCTION 


When flames are stabilized on bluff bodies, the interaction of a cool 


combustible mixture with its hot combustion products is of great importance 
for the stabilization mechanism (Zukoski & Marble 1955; Cheng & Kovitz 
1958). An idealization of the problem was first carried out by Marble & 
\damson (1954). ‘They consider a non-viscous, perfectly insulating, 
semi-infinite partition separating a cool combustible -mixture from its hot 
combustion products. At the trailing edge of the partition velocity, 
temperature and concentration of combustible are uniform in each half- 
plane. ‘Temperature distribution and combustible concentration are, 
initially (ie. at the trailing edge), step functions while the velocity 
distribution may or may not be a step function. Their analysis then gives 


the temperature, combustible and velocity distributions in the neighbourhood 
of the trailing edge of the semi-infinite partition. Of particular interest 
is the distance downstream from the trailing edge at which a local 
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temperature maximum, due to heat release by chemical reaction, first 
appears. 

The present analysis is concerned with the identical problem except 
for one important difference. ‘The partition between the two streams is 
now taken to be of finite length and viscosity is included. As a result, 
the initial velocity distribution is no longer uniform but of the Blasius type. 
It will be seen that this significantly alters the distribution of stream properties 
in the immediate neighbourhood of the trailing edge as compared with 
Marble & Adamson’s result. 


2. ASSUMPTIONS CONCERNING THE MODEL AND ANALYTICAL FORMULATION 
OF THE PROBLEM 

In this paper a detailed discussion of the assumptions on which the 

model is based will be omitted. They have been fully discussed by Adamson 

(1954). In essence, the equations to be presented for solution describe 

the following flow, which is schematically shown in figure 1. In the upper 
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Figure 1. Schema of model with coordinate system. 


half-plane there is, initially, a cool combustible with temperature T,, 
density p,, and free stream velocity u,;. ‘The lower half-plane consists, 
initially, of fluid which is chemically inert with respect to the upper stream 
(say its combustion products). Its temperature is 7, > 7’, its free stream 
velocity is u,;,; = u,, and its density is p,;,; = p; T,/T;, since the pressure is 
assumed constant throughout the field. ‘The specific heat C,, and molecular 
weight VM of each component (combustible and combustion product) are 
equal, with C,, taken as constant. The transport properties, dynamic 
viscosity jy, coefficient of thermal conductivity A, and binary diffusion 
coefficient D of each component are equal and vary as though the molecules 
were Maxwellian with the Eucken correction for A holding. From the 
above it follows that the Prandtl number Pr = C,, A/u and Schmidt number 


F.M. E 
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Sc = p/pD are constant throughout the field. Also, pu = constant. Heat 
release by chemical reaction is assumed to be first order, that is, the heat 
release O per unit volume per unit time is assumed to be given by the 
equation 
QO = Be oie ART 

where AH = heat release per unit mass of combustible, p = local density, 
mass per unit volume, A = relative mass concentration of combustible, 
mass of combustible per unit mass of mixture, 7 = characteristic chemical 
time constant, R = universal gas constant, A = activation energy, a constant. 
Finally, the flow is assumed steady with a laminar mixing region in which the 
usual boundary layer approximations hold (Goldstein 1938). 

The governing equations are obtained from conservation of mass, 
conservation of momentum, conservation of energy, and conservation of 
combustible. Since pz = constant the Howarth transformation (Howarth 
1948) may be used to uncouple the continuity and momentum equations 
from the energy and combustible equations. ‘The equations to be solved 
can then be written as 


u,+v, = V, 
uu,+VU, Uys 
ub, +o = - oe 4+. BKe-%/6 
(2) 
uK,+vK, = —K,,—CKe-™ | 
Sc 
where B= 4/AH|u,yC,7,,, C =4l/u,7. 0= T/T, 8, = A/RT,,, and 


lis the length of the flat plate. ‘The symbols wand 7 denote non-dimensional 
axial and transverse velocity components in the incompressible plane 
while x, y are the corresponding distance coordinates, with origin at the 
trailing edge of the flat plate. Their definitions are taken from Goldstein 
(1930). 

Equations (1) have been solved by Goldstein (1930) for the case of flow 
in the wake of a flat plate with symmetric Blasius velocity profiles at the 
trailing edge. ‘This solution will be used in solving (2). 

Goldstein’s solution for u, v is expressed in terms of the independent 
variables 


I 


& y! 3. y y 3§ 
These variables will also be used in the solution for 6and AK. Furthermore, 
since the Blasius initial condition is in terms of the Blasius series uo(1), 
valid only tor y > 0, Goldstein’s solution for u,v is valid only for n > 0. 
It is therefore necessary to construct solutions for 6 and K in the two 
half-planes separately and join them smoothly at y = 0. (Quantities in 
the upper half-plane will be denoted by an overbar.) 

Another consequence of the Blasius series initial condition is that 


Goldstein’s solution is expressed in two forms, one valid for small y, 
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the other valid for large y. Therefore, solutions for 6 and K must also be 
constructed in the same way. 
The boundary conditions for 6 and K are 
lim, 9(£, 7) = lim, 0(€, 7) = T,/T,, = A, 
= lim, O(€,) = 1, 


(€,) 
lim, K(é, 7) = lim, K(£, 4) = 1, 
(€,7) 


lim, 6(€, 7 


lim, K(£,7) = lim, A(&, 7) = 9, 
where 
lim, = §-0, 74>, y=const., 
lim, = 7>®, yo, 
lim, = €>0, 1->%, y=const., 
lim, = 7>®, y>o. 


3. SOLUTION FOR @ AND AK FOR SMALL y 


For small y, solutions of the form 


x 

4,0) = 2 A,(n)é", 
pry (3) 
x 

K(é,9) = X K,,()&"s 

are assumed, with similar forms in the upper half-plane. 

Since the initial condition as €--0 with y = constant is independent 
of y it can be shown that the first term of each series must be independent 
of € and all terms with nm > 1 must vanish as £0, y = constant. In terms 
of the assumed solutions (3) this implies 


lim 6, = A, lim 6, = 1, 
i i (4) 
lim A, = 1, lim K, = 0, 
and i = 
lim 6,,/7” = lim @,/n" = lim K,,/4”" = lim K,,/n" = 0 (5) 
for n > 1. 


For large y the form of solution must be such that at some intermediate 
value of y the two solutions join smoothly. Its form will be discussed after 
the solution for small y has been obtained. 

Goldstein (1930) gives, for small y, 


3 3/+1 ’ 
itis (6) 
v= 1 > C3) E3k—1 


where 


E2 
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The functions fy, fs, f, and their derivatives were calculated by Goldstein 
and have been recalculated for smaller mesh size in the present work. 
Substituting (3) and (6) into (2), and equating the coefficients of like powers 
of € to zero in the usual manner, we obtain a series of equations for 6, (7) 
and K, (7): 

0° T Ag fy 9; S Af, 9, Pri 


K" +yofoK, —y.f, K, = ScH,, ) 
where 
\ wy, ¥ SSC, fe ©, 1, Zc) 
L.-i,=, = 8, «9%, 
I, - 9BK, e~ 4’, 


and for $2 3 
i 9B > D;K,+4 > ((¢-3)f5,%-s— (31+ 2) fa, 9-3} 


with H, identically the same except for replacement of B by —C, 4-5 
by K,, and 6, by A, ,. The coefficients D; are obtained from the 
expansion 

ee . 6.6 
eW4 — 5 D (n)& = e@77aio + ok FF 

i=0 01) 

Similar results hold for the upper half-plane. 
By the method of variation of parameters a formal solution of (7) may 


be written as 


A.(n) = a, f(y) +6, fo'(n)—Pr | L,(m)expsa, | fol) ni | c 

. é \ i ~ . * r : : 

< f(a) fP () fi? @fz'@)5 4, (8) 

where a,, 6, are arbitrary constants and f\, f5" are linearly independent 

solutions of the reduced (homogeneous) equation. A similar result holds 
for K (7). 


For small and large », fy behaves as 





; 
P ) a2 |] 923 7 
\ _ — 2B 


3! 5 i le’ (9) 
hy = 4 (1) y d9)° T Ol(y u do) Pexp| i 5% (7 t do)*t], 
respectively, where 8, = 3:67869, x, = 5-97708, 55 = 0-3408 (Goldstein 

1930). 
For small » a solution of the reduced form of (7) is given by 
OY) = a(l+ySyPtylntt ...)+b,(y + vonF + yon? + ...) 
= a, fi +6, fs”, (10) 
where the y\*) are known functions of 89, A,, and A,. ‘The meaning of a, 
and b, in (8) is now clear. a, is the value of 6, at 7 = 0, while 4, is the value 


of @, at ny = 0. 
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4. SOLUTION OF THE REDUCED EQUATION FOR LARGE 7 
The boundary conditions (4) and (5) are all applied as 7 > 2%. Therefore, 
a, and b, cannot be evaluated from the general solution (8) until the 
asymptotic form of this solution is known. Since /f\*) and f$*) cannot be 
expressed analytically for all » the asymptotic form of (8) cannot be 
obtained directly. 
Using the asymptotic form of f) and neglecting the exponentially 
small term, we find that the equation (7) for 6” becomes 
A" + aa?G” — sac8 = 0 (11) 
where a = a, Prand o = 7+. The solution of (11) gives 6\” as a function 
of o. It can be shown that only exponentially small terms will be added 
to &\ as 7 x if the exponentially small term in f, is retained. 
The substitutions 


3 13 2 13 2 
IW IW : ef J@ 
o=(—) , OM =(—) eH {— 
a a a 


transform (11) into 
2 F f 1 Y 1 1 
d Hi Se l 3 ( l ) 4 =i6 , W. 3 0), 


— — + + 








(12) 


du)” | + WwW an 
which is a form of Whittaker’s equation (Whittaker & Watson 1952), 
In terms of #6 and o the asymptotic form of its solution may be written 








ey 
+ B,o8{ 14 
\ 
2 F tele bs— 34+ its (“bs 5+ tes (5+ bt 7) 
a) n!(—}ao%)" J 
(13) 


where A, and B, are arbitrary constants. A similar result holds for the 
upper half-plane and for K‘”. 


5. DETERMINATION OF ARBITRARY CONSTANTS FROM BOUNDARY CONDITIONS 

Although a particular solution to the asymptotic form of (7) has not, as 
yet, been obtained it can be seen from (13) that A, will not be fixed by the 
boundary condition for the lower half-plane alone. However, B, must 
be determined since the second term of (13) becomes of order o* as o> ~. 
It will be seen that as o— « the particular solution of (7) becomes 
exponentially small in the lower half-plane and can be neglected in the 
upper half-plane due to the assumption of zero heat release by chemical 
reaction in the cool combustible. This, therefore, implies that for s > 1 
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s) 


Six constants remain to be determined. ‘They are a,, b,, A,, d,, b., 
and A,. ‘Their determination stems from the matching of value and slope 
of 6, and @, at two values of 7, namely 7 = 7 = 0 and 7 = 7 = 9, 7, being 
some intermediate value of 7 and 7 where both solutions of 6, and 6. are 
assumed valid. 


The matching of value and slope at 7 = 7 = 0 immediately gives 
I ] } : 


a, = Ay, 
- (14) 
b, = —b.. 

The matching of value and slope at 7 = 7 = »,, gives 

a, D,,+6, Do, + Dy, = A, Dy, + Des, 
a, D,, +5, D3, + D;, = A, D+ D;,» P 
: ; - % (15) 

a, D,,—6, Do, + Dog = Ag Dy.+ Dogs 


a, D,, —b, D3, + Dj, = A, D,,, + Dt, 


8 


where the D, and D! are functions of 7 evaluated at 7 = 7 = »,,. Thus six 
equations are available for the determination of the six arbitrary constants. 


6. "THE FIRST FIVE TERMS OF THE SOLUTION FOR SMALL ¥ 
(i) For large 7, from equation (13), 
6, ~ Ayo exp(— 4a0*)< 14 > a%o-* \ +1, (16) 
where a” is a constant given by (13) and B, = 1 because of the boundary 


condition (4) and, for large 4, 


6, ~ Ayo- 2 exp(— $a0*)d14+ ¥ ao" $4 A, (17) 


l 


since B, = A. 
Carrying out the matching technique described in the previous section, 
it can be shown that 
| s(1 + A), by (1 \) [Day D5, Dp BD.) Ay = Ay, (18) 
where the D symbols are as defined by (15). 
Since the D’s do not depend upon A it can further be shown that 
1+A 1—A 1 + AY] 
: ‘i es * 
G4(; Pr A) = “se + Toa 6,(y, Pr, A ) D - | (19) 
where A* is a particular value of A. ‘The solution for Ay is obtained in 
the same way with AK, = } when 7 = 0. 





The function 6 (y,Pr,A*) has been integrated numerically for 
Pr = 0-91 and 1-00 with A* = 0-286. Numerical integration is discussed 
in $9. Figure 2 shows 6, and K, for Sc = 1-00. 

Note that 


Ko (1 — 6%) (1 om FP}, (20) 
(ii) In solving for 6, through the matching technique it can be shown 
that a non-trivial solution for a, and b, will be obtained only if the Wronskian 
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FPFO” —fIP’fS? vanishes. Since f}!’ and f}’ are linearly independent solutions 
for 6, this implies that the only solution for 6, satisfying the boundary 
conditions is that for which a, = 6, = 0 or 0,=0. The same result holds 
for A,, namely, A,= 0. 

(iii) The equation for 4, is 


05 + 2f) Pr 6, —2f; Pré, = —9BPrK,e ~~, (21) 
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Figure 2. Solution for 6, and Ko. 


4, is the first factor in the solution for 6 which contains the effect of chemical 
reaction. This is seen from the form of J,.. With chemical reaction absent, 
B=0(0. Then #,=0, by the same argument as was used to show that 
6,=0. By direct substitution it can be shown that if @,(n,...,B*) is a 
satisfactory solution of (21) for B = B*, then, for any B, 


8,(7, Pr, Sc, 8,,.A, B) = a als Pr, Sc, 6,, A, B*). (22) 


‘Thus the dependence on the parameter B has been extracted. Furthermore, 
using the matching technique of $5, with exp(—4@,/A) <1, it can be 
shown that 
@,0c e7 (23) 
for values of » where @, = 1. 
It is interesting to note that if Sc = Pr, 


K.(n, Pr, 8,,A,C) = — = 6.(n, Pr, Sc, 0,, A, B*). (24) 
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K, is, physically, the first term in the concentration solution that gives 
the effect of loss of combustible due to chemical reaction. 

6, has been obtained by numerical integration for B = 2-58 x 104, 
6,, = 23-96, Pr = 1-00, and A = 0-286. It is shown in figure 3. Outline 
of the numerical method is deferred to § 9. 
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Figure 3. ‘ Temperature distribution ’ 4, due to heat release J, 


From the dependence of /, on Ag, it can be seen that /, becomes 
exponentially small as 7-» «. ‘Therefore, a particular solution of (21) 
for large » would also be of exponential order. In the upper half-plane 
6, A, so that J, is of the order of Bexp(—4,/A) < 1, since chemical 
reaction is assumed to be negligible in the cool combustible. ‘Therefore, 
a particular solution of (21) for large 7 would be of order Bexp(—4@,/A). 
These observations show that for the boundary conditions to hold, 


B, = B,=0. Thus, for large y and a = 4, 


x Ip, 
0, ~| Aro 1+ Sa2o-an| 4 FBP TE o—rapg—3_ 1H 9-6 4 3 | ta 
L a? } a 
(25) 

(iv) Since #, = A,= 0, it can be seen from the definition of J; that 6, 

is really a term due to mixing without chemical reaction. Furthermore, 
I, and 1, are antisymmetric about 7 = 7 = 0, because 4’, and @’ are anti- 
symmetric about 7 = 7 = 0. This implies the antisymmetry of @, and @, 
about 7 = 7 = 0. ‘This conclusion may also be seen from the matching 


technique of §5. 
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Finally, because of (19) and the dependence of J; on 4), 
1A 65, Pr, A* (26) 
a= 0,(1, Pr, AY), 2 
63 is exhibited in figure 4 for Pr = 1-00 and A* = 0-286. 
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Figure 5. Solution for 6,, showing J. 


(v) 6, is the last term computed with chemical reaction. A particular 
solution of the asymptotic form of the governing equation is exponentially 
small since J, becomes exponentially small when 7 ©. ‘Therefore, as 
before, B, = 0. The same argument as used to show that B, = 0 will 


show that B, = 0. 6, is exhibited in figure 5. 
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(vi) The solution of 6(¢, 7) for large » and small y will now be collected 
since it will be useful in the discussion of the next section. Ifthe asymptotic 
forms of 6,(y) given by (13), together with the asymptotic form of 
particular solution to @,, be substituted into (3), the asymptotic form of 
H(€,n) is obtained. Using #, and @, only, with Pr = Sc, i.e. a = b, and 


G =17T 9%, 
A(E, n)~ 14 [Ao{y “2 2657) 34 367 7) +...) 
+ E27A,{n—4 — 46,7-° +...) + ... Jexp} - La(y? + 385 7? + 387 y +63)! (27) 


after making use of B. = B.=0. A similar result holds in the upper 
half-plane and for K(€, 7). 


7. SOLUTION FOR LARGE y AND ITS MATCHING WITH SOLUTION FOR 
SMALL ¥Y 

It was pointed out in § 2 that Goldstein’s (1930) solution for uw, v is split 
into two forms, one valid for small y, the other valid for large y. Accordingly, 
the solution for @ and K must also be so split and made to match at some 
intermediate y which is really the outer limit of validity for the small-y 
solution. 

The form of solution for small y at its outer limit of validity is obtained 
by making 7 large. This form, for 6(€,7), is given by (27). ‘The new 
solution valid for large y will be forced to match this form when y is 
reasonably small. 

For large y it will be assumed that the terms measuring the rate of heat 
release due to chemical reaction and the rate of consumption of combustible 
are small compared with convection, conduction, and diffusion terms. 
This assumption implies that the region in which Goldstein’s wake solution 
for large y applies is one where pure mixing, without chemical reaction, is 
predominant. The governing equations (2) are then simplified to 


es > 1 s ok . Y.—lK” 2 
u0.+v0, = Pr-6,,, uk ,+vK, = Se"K,,, (28) 
with identical equations in the upper half-plane. 
Goldstein’s velocity solution, for large y, to be inserted here, is 
hs. us. 
/ Cp! 4 6812 2 ER ES 
u Py + op, +E 71 $ 3! rane 
ail (29) 
us us. 
- L¢ 2b 2¢ 172 3 ¢0 73 _ 
. a , 35 C ’ 
|, 
yp? a , us Z. : ” tbe Ps 3 >) 
where hy = p(y), y, = Au,, xy = 5740 71 Ty ln — ilo» 


and A = 36, with u(y) representing the Blasius velocity profile at x = € = 0. 
For given sufficiently large y, € can be made arbitrarily small. In so 
doing, v is made arbitrarily small, since 7 = y/3€. Considering only terms 


independent of €, equation (27) becomes, as €>0, 


A(E,n) ~ 1+ Ag(q*— 28) 7-2 +... )exp{ — Ja(n? + 38, 1? +352 +83)}. (30) 
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The solution of (28) using the u, v given by (29) must reduce to (30) 
for €--0 and sufficiently small y. This suggests assuming a solution of 
the form 

A(f,9) = 1+N > O,(y)é"*? exp{O,(y)E3 + O,(y)E* + O.(y)EM (31) 
in the lower half-plane where ©,, 0,, 0, and ©,(y), with r = 0, 1, 2, ... 
are, as yet, undetermined functions of y, and N is an, as yet, arbitrary 
constant. Substituting this assumed solution and (29) into (28), and 
equating the coefficients of like powers of é"exp{O,,€-3+...} to zero, one 
gets a series of equations for 0,. 

The first four equations are 











1 : ) 
x 2440, = 0, | 
r | 
se gil 1 
Pp, = - O, +, —(+-u4, = 0, 
r J Ja a ») | 

a 102 2,,0 hs, © i. & De 

O+- 0+ — +e ae + pa + oh +e = Oe 

Pr °' 30! BG, 3° "216 * 3° 32 [ 
2 ? 3 ? | 
— (log @,)’ +4 = d .% ae ee oy _2 Be = wy | 4 | 
> QO, dyi3** * 32 | | 


Although the equations after the first one remain first order and linear, 
they rapidly become lengthy. 
Solving for ©, one gets 
1 | 2 
O,= — —<G,—-3Pr | fi? dy} 33 
By") (33) 


where G, is an arbitrary constant. From the definition of ¢/, 


wb, = u(y) =a, ytay'+..., 


) 


where a, = ta, a, = — }a?/4), 
and a = 1-32824. Goldstein defines a,,,., = 3***a,,.,, for r = 0, 1, 2, .... 
hen @, = a,/9 = a/9Pr, 


since a = x, Pr from (11). Therefore, for sufficiently small y, 


OQ, ~ —Pr-4G,— Pr 2al2y32)2, 


a 


Pr af y\3 
O,~ —-—a ea = 
ie 3\3 


This is exactly the form required for the matching. ‘The required solution 


lf G, = G, 


of ©,, is then 
esd bP | Win dy | (34) 
The solution for ©, can be shown to be 


P - 23 
©, = —3Pr Ay) . bi : dy i o, | ih) . dy | ’ 
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where G, is an arbitrary constant. Setting G, = V0, the behaviour of 0, 
for small y is as required, namely, 
©, ~ —ado(y/3)?. 


‘The required solution for ©, is then 


ui 
~—" 


©, = — 3Prd,y,1? | WV? dy, (3 


7 


since A = 36,. 
The required solution for ©, is 


(-) -2 Prai| yb U2 | bi? dy t+ Wi | (36) 
For sufficiently small y it becomes 


©). ~ —abo(y/3) 


as required for perfect matching with (30). 
The solution for ©, for all y is lengthy but its behaviour at small y can 
be checked to be 


as required for the matching. 
Finally, the constant V must be 
N = 9A,exp(— 463). 
Terms of higher order than ©, were not calculated. 

Thus a solution of @ for large y has been obtained and matches the 
solution for small y in the following sense. For a given sufficiently large 7, 
€ is made so small that only terms independent of € are important in the 
asymptotic solution stemming from the region of small y but large 7. ‘The 
solution for large y matches this form as y gets sufficiently small. 

The series (31) does not permit matching with the series solution for 
small y to terms involving €? and higher as appeared in (27) for large 7. 
Attempts at modifying the formal expansion (31) to obtain a formal matching 
to higher order powers of € have so far been unsuccessful. ‘This is in sharp 
contradiction with Goldstein’s velocity solution where perfect matching 
is indicated. 

It is rather fortunate, as shown by the computed results, that the physical 
problem of present interest is in the region of very small € and a value of 7 
such that the solution for small y applies. Therefore, the matching 
procedure, unsatisfactory as just described, recedes in importance except 
for mathematical interest. 


8. LOCATION OF FIRST LOCAL TEMPERATURE MAXIMUM 
For a pure mixing interaction between hot and cold flow the temperature 
distribution in the wake follows a uniformization process. No local maxima 
can appear since there are no heat sources in the flow field. With heat 


release by chemical reaction this is no longer true. The distributed chemical 
reaction raises the temperature throughout the field and with the assumed 
rate law given in $2 a local maximum in reaction rate occurs in the hot 
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stream thus causing a local temperature maximum to appear there also. 
This local maximum then rapidly propagates upward into the cool 
combustible. 

With the temperature distribution given by 6(£,), the locus of 
temperature maxima, (7), is given by 6, = 0. 

Prior to formation of the initial local maximum, 6, <0 for E< &, 
where &; is the location of the initial maximum. After the ‘maximum forms 
there must be a value of 7 at which an inflection occurs, i.e. #,, = 0. Above 
this point #,, <0, and below this point 6,,>0. As §—€—+0 the 
inflection point coalesces with the point where 0,=0. Therefore, the 
initial maximum is given by the conditions 


se nf Nie 
The numerical results indicate that it is sufficiently accurate to use the 
first two terms of 4(€,7) for calculation of €;. From (37) 


6, /0, = 6, 6; at 7 = n,, 


and 
é&=(-0 gi)! - si (38) 
From (19), (22) and (23) 
4’ (n, A) 1-A 6'(n, A*), 93(y, 8, B) = B exp(4* — 0,,)03(n, 0*, B*), 
1—A* " 2 - B* ce, 7 


where it has been assumed that the initial maximum occurs where 4, = 1. 
Putting these results into (38) one gets €; in terms of a known €*. With 


the restrictions that Pr = Pr*, Sc = Sc*, this gives 
1—A B* ag |? A 
go fie’ Bae, & (39) 
1—A* B 


This simple formula (39) gives explicitly the dependence of the critical 
distance €, on all the important dimensionless quantities, namely, the 
temperature ratio of the two streams A, the activation energy 4, = A/RT,,, 
and most interesting of all, the parameter B = 4/AHj,/u,7¢, 7,, which 
combines the fluid mechanical properties /, u, with the chemical kinetic 
properties AH,,/c,, 7, and 7+. 

‘The effect of each individual physical quantity can be obtained easily 
from this formula. A detailed discussion and the implication concerning 
the appearance of a flame has been given in a previous paper (Cheng & 
Kovitz 1956). 

The following important conclusions may be mentioned. 

(i) ‘The presence of a viscous wake shortens, by an order of magnitude, 
the distance €; where the temperature maximum occurs, as compared 
with the result without the initial wake (Marble & Adamson 1954). 

(ii) The apparent fluid-mechanical variable in determining the distance 

€, is the shear cu/dy at the trailing edge of the plate. 
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iii) The variation of €, due to thermo-chemical-kinetic properties is 

(111) g, prop 
overwhelming compared with the variation produced by the 
fluid-mechanical properties within practical ranges. 


9. NUMERICAL METHOD 

Since Goldstein’s velocity solution involves a numerical tabulation 
of the functions f,, fy, f, it is, of course, necessary to determine 4), 4,, 42, 
numerically also. ‘The system of equations for the 4, being linear, it has 
been possible to extract the parameters A from 4) and 43, and B from 64. 
That is to say, a tabulation of 6, and 4, for a given A yields 6, and 4, for 
all A from (19) and (26), while a tabulation of 6, for given B yields 6, for 
all B from (22). In $10, it will be seen that when chemical reaction is 
absent A may be extracted from all the non-vanishing coefficients 
6,,, r= 0,1,2,.... The method of trial and error is preferred to the 
superposition of solutions to meet the boundary conditions at two different 
points in the numerical integration for higher accuracy. 

The numerical integration of 4 and Ky is straightforward since the 
initial values at 7 = 7 = 0 are specified as in §5. The initial value of 6, 
is also fixed from the result of $6. These functions were obtained from a 
single trial and error procedure. 

First approximations to 6, and 65 at 7 = 7 = 0 can be obtained from 
the matching technique of §5, but their precise determination was made 
from a double trial and error procedure. ‘The single pair of values of 
6,(0) and 65(0) which give correct behaviour of #, and @, were rather quickly 
localized using a graphical interpolation method. Details are given by 
Kovitz (1956). 


10. SOLUTION WHEN CHEMICAL REACTION IS ABSENT 
The absence of chemical reaction implies B= C=0. Within the 
approximations used, this uncouples the dependence of 8 on K as seen 
from (28). From the general results in $3 the inhomogeneous term then 
becomes, for s > 3, 


I, a At 3 iy 1, 3 - (3/ t 2) fa 9, fe (40) 


It has been shown that 6,=0. Therefore, /,=0. By the same 
reasoning as was used to show that @,= 0, it can be seen that 6,= 0. 

With B = 0, I, = 0, as seen from the definition of J, in $3. Therefore, 
as previously shown for 6,, 4,= 0, which implies /;= 0, which in turn 
implies 6; = 0. 

From equation (40), 
[,= — 5f,0/ ae fi, 6. — She 6’. 


From Goldstein (1930) f, =f, = 0 at » =0. Furthermore, it has been 
shown that @,(7 = 0) = 0 and that @’, 6,, 4, are geometrically antisymmetric 
about »=7=0. Therefore, J, is antisymmetric about 7 = 7 = 0. 
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A single trial and error numerical integration yields @, at » = 0 such that 
the boundary condition 
lim 6,/n* = 0 


n> 


is met. 46, is shown in figure 6 for A = 0-286, Pr = 1-00. 
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Figure 6. Solution for 4. 


From the form of J, and equations (19) and (26) it is clear that 
1—A 

~ A* 

From (40) it is seen by induction that 7, depends only on (s—3)6,_, 

and 6’, for s >3. Therefore, only J;, (r = 1, 2, 3, ...) is non-vanishing, 
which implies that only 43, (r = 1, 2, 3, ...) is non-vanishing. The general 
solution for all A in terms of a solution for A* and given Pr is then 

; > 1 1—A f 1 *) 1A dy Ae 
6(€,7, Pr, A) = 3(] + A)4 1-A*' -S(1+A )+ o(7, Pr, A )+ 

+ €36,(n, Pr, A*) + €90,(n, Pr, A*)+...}. (42) 


6,(y, Pr, A*). (41) 





6m, Pr, A} = 


11. CONCLUDING REMARKS 

The region of convergence of the series solutions cannot be specified 
analytically. However, the numerical results do show that with chemical 
reaction the magnitude of the coefficients @,,(7) increases rapidly. Since 
this is essentially a linearized solution only, the region where the effects 
of chemical reaction are small compared with pure mixing effects can be 
considered. This is delineated by 0 < & < &; where &€; is again the value 
of € at the appearance of the first local temperature maximum. A detailed 
computation (Kovitz 1956) using the first four terms of 8, with azomethane 
as combustible, shows good convergence for € = €, = 0-048. Considering 
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only the dependence of €; on B, given by (40), it can be seen that if the 
series converges for &*, it will converge for €,, since €/ 6, is approximately 
independent of B. 

If chemical reaction is absent the coefficients 6, grow less rapidly. 
‘The magnitude of /,, and so the magnitude of 6,, depends upon the magnitude 
of Goldstein’s f,,. It may thus be estimated that the convergence without 
chemical reaction is of the same character as the convergence of Goldstein’s 


solution for u, v. 


‘The authors wish to acknowledge the invaluable aid of Mr R. A. Goerss 
in setting up and performing the computations. ‘This research is part of 
the work under OOR DA-36, ORD-2183, sponsored by the Office of 


Ordinance Research, U.S. Army. 
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SUMMARY 


A numerical solution has been obtained for the starting flow of 
a viscous fluid past a circular cylinder at Reynolds numbers 40 
and 100. The method used is the step-by-step forward integration 
in time of Helmholtz’s vorticity equation. ‘The advantage of 
working with the vorticity is that calculations can be confined to 
the region of non-zero vorticity near the cylinder. 

The general features of the flow, including the formation of 
the eddies attached to the rear of the cylinder, have been determined, 
and the drag has been calculated. At R = 40 the drag on the 
cylinder decreases with time to a value very near that for the 


steady flow. 


1. INTRODUCTION 

The general features of the flow of a viscous fluid past a circular cylinder 
are known from experiments. At low Reynolds numbers a_ steady 
symmetrical flow exists. At Reynolds numbers above about 40 the 
symmetrical flow is unstable and periodic oscillations known as the 
Karman vortex sheet appear. Numerical solutions of the steady flow 
have been obtained at Reynolds numbers 10 and 20 by Thom (1933) and 
at R = 40 by Kawaguti (1953). 

In this paper the two-dimensional flow of a viscous fluid started 
impulsively from rest and moving perpendicular to the axis of a circular 
cylinder is investigated. Apart from its intrinsic interest, this flow is 
expected to throw light on the high speed flow around yawed bodies of 
revolution. The method employed, which is the integration of Helmholtz’s 
vorticity equation, is that which was used by the author (Payne 1956) for 
the starting and perturbation of a two-dimensional jet. The preliminary 
results, described below, were obtained for the starting flow past a circular 
cylinder at Reynolds numbers 40 and 100. In these calculations no 
asymmetry was introduced. With the aid of a more powerful electronic 
computer, it is intended to proceed to higher Reynolds numbers and also 
to investigate the unsymmetrical flow. 


F.M. F 
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2. EQUATIONS OF MOTION AND BOUNDARY CONDITIONS 


f Helmholtz’s vorticity equation with allowance for 


Use is made 


viscosity, 


Cw al uw) o(vw) 2 O*w fs *W) 
: oe a T eens ae a. ae a T ag Fa ( 1 ) 
et Cx cy R |ox*® = dy? 


R being the Reynolds number, to find the vorticity w at successive small 
intervals of time. (‘The cylinder has unit radius and the velocity at infinity 
is unity.) ‘The velocity (u,v) is found by adding to the velocity of the 
classical inviscid steady flow the velocity due to the vorticity, which is given 


by 





Lt an 
u(x, y, t) i | | aX, Y,2) 7 X¥F+(y 2 hailed 


l ¢ (x —X) 
u(x, y,t = wf X, Y,t <5 
ssi i dia A Gy X)?+(y—Y) 
The boundary condition that on the cylinder there is no flow perpendicular 


to the surface is satished by introducing image vorticity, while the no-slip 
condition causes vorticity to be generated on the cylinder. ‘The vorticity 


, dXaY. 





is ditfused and convected away, giving a region of non-zero vorticity near 
the cylinder. ‘The calculations are confined to this region. 

It is convenient to use coordinates € = logr, where r is the radial distance, 
and 6, the azimuthal angle. Helmholtz’s vorticity equation becomes 


a - ~ >) 2 A9 


( ( C°W O°W 17 
—F (rl w)-+ ay (rI w) = 
cc ( 


apy Tt ApB > (3) 


where U and F are the radial and transverse components of velocity. 


3. FINITE DIFFERENCES 
The plane is covered by a mesh of points c RA, 0 ja, where k and i] 
are integers and Ais aconstant. From equation (3), and replacing the space 
derivatives by central difference formulae and using a forward difference 
for the time derivative, a first approximation to the vorticity at the next 


step in time is 


Q = wr + At F(UZ,, Ver,, wey) 
where 
F(l V ) [ l 
of | Ix J» ] ) j - n ‘ Ww r; w 
KJ ia 2r? A ; ; 
Br ee —r,I We iaitt 
z 
+ = Ww, «a +y, + W) tw! — 4a}. 
Rr? A? 1,3 LJ , we 





The velocity at time (7+1)At is then evaluated from the integrals (2), as 


approximated by equations (6), below giving 














~ 
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Then a second approximation to the vorticity at time (m+ 1)At is 


(n l nt l n) A Tin 7(n+1 n+ly\ 
Wy. 5 = FQ y 4 4 {wv ‘Te 2 At i (l K,J 1) I KJ ON )}. (4) 


3 
Since the vorticity at each lattice-point is A~* times the circulation around 
a quadrilateral of area A* centred on the point, and this circulation is VA 
if the transverse velocity is V just away from the surface and 0 on the surface, 
the vorticity at points on the cylinder is given by 








afte? = VEYA. (5) 
The velocity due to the vorticity is obtained from the formulae 
s, Ver) _— 2 
ae = > > (A, m,j—l A, rim iD" n 9) =f *, 
A oe l =a A ; ' : ; i 
F (6) 
‘eer ~ } RB ay, Om+n 
RT 2k igh Birmsrm Brrmiat Fe Jin Ons | 
where 
A sin bA , 
Aj,,=-7> if a~0 or bX 0, 
-” 2m exp(aA) + exp(—aA)—2cosbA # 
A exp(aA)—cosbA o 
eae a ifa+0 or b¥<0, 
”  2nrexp (a\)+exp(—a@A)—2cosbA 
and Ago = Boo = 0. 


4. DRAG 


The total drag on the cylinder is equal to the rate of decrease of momentum 
of the fluid, from which it follows (Phillips 1956) that 


d 7, 
D=-- 5. || 


yw dxdy >. 
Chae 


ry*2> J 
The drag coefficient is Cp = D/(}pU%,2r)) (where rp = U, =1 in our 
calculation, ry being the radius of the cylinder). Using the coordinates € 
and 6, and replacing the double integral by a double sum, the drag coefficient 
at time (n+ 4)At is found as 


1 at Wh n+l 
Cza= = — SN *% r? sin Nas) A2 ; 
D At BS o .* JAW; 


For Reynolds number 40 the drag coefficient thus calculated was 
initially 3-00, decreasing rapidly to 1-91 after unit time (the time in which 
the fluid at infinity moves a distance equal to the radius of the cylinder) 
and thereafter decreasing more slowly (see figure 1). At ¢ = 6 it is about 
2°, above the value C, = 1:6177 calculated by Kawaguti for the steady 
flow at R = 40. For Reynolds number 100, the drag coefficient was initially 
1-20, decreasing to a minimum of 1-00 just before t = 1 and rising to 1-10 
at ¢ = 3. It is easily shown that for all Reynolds numbers there is also an 
impulsive drag coefficient of 7 when the flow is started. 

F2 
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For small t, the boundary layer is very thin so the value of Cp given by 
equation (8) will not be accurate. The coarse mesh will not satisfactorily 
cover the flow within the boundary layer itself. In this range of ¢ it is 
necessary to use the boundary layer solution of the equations of motion. 
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Figure 1. The variation of the drag coefficient with time for the starting flow past a 
circular cylinder. ‘The value obtained by Kawaguti for the steady flow at 
Reynolds number 40 is also shown, together with the analytical solution for 
small t. ‘The boundary layer solution is not valid beyond the instant at 
which separation occurs; however, beyond separation it is shown as a broken 
line suggesting a method of joining it to the numerical solution. 


To the first approximation as t-> 0 (Goldstein & Rosenhead 1936), 


V —2sin0 . | exp(— x") dx, where 7 cent 
—_ ee 
V (arty “PSTD 
From equation (7) it then follows that 
at ° aris | 
i a j | of | ; wr* sin 6 nad | 
- vs (27 + 84/(zvt) + O(t)} 
dt ' 
~4,\/(27/Rt) as t-0. 


It is interesting to note that skin friction and 
pressure contribute equally to the total drag for small t, the respective drag 
being 2, (27/Rt). Since the boundary layer ultimately 


This is also shown in figure 1. 


coefficients 


separates, this value of the drag will not be reliable beyond (at most) 
t = 0-32, when separation occurs (according to the boundary layer solution). 
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Unsteady viscous flow past a circular cylinder 








Figure 2. The symmetrical starting flow past a circular cylinder at Reynolds number 
100 at time 2. In the upper half is shown the vorticity distribution and in the 
lower half the streamlines, the direction of flow being from left to right. The 
vorticity at the rear of the cylinder causes a reversed flow near the rear stanga- 


tion point. 
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Figure 3. The symmetrical starting flow at Reynolds number 100 at time 6. 
two eddies attached to the rear of the cylinder have increased in size. 





The 





Figure 4. The symmetrical starting flow at Reynolds number 40 at time 6. 
vorticity spreads further out laterally at this lower Reynolds number and there 


4 


is a greater reduction in the velocity near the cylinder. 
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5. RESULTS 

The symmetrical flow past a cylinder starting impulsively from rest 
has been calculated for Reynolds numbers 100 (figures 2 and 3) and 40 
(figure 4). ‘The space mesh size used was A = ed = 0:20944 so that 
there were 30 mesh points round the cylinder at intervals of 12 degrees. 
‘The time interval was taken as At = 0-1. 

The vorticity generated at the cylinder is transported round towards 
the rear stagnation point. ‘This vorticity induces a reversed flow on the 
cylinder, the velocity component at the mesh points nearest the rear stagnation 
point changing sign at ¢ = 1-5. The two eddies attached to the rear of the 
cylinder appear and increase in size. At ¢t = 6-0 they have spread 1-5 radii 
downstream (at both Reynolds numbers 100 and 40), which is about half 
the length of the standing eddies in Kawaguti’s steady flow at Reynolds 
number 40. 

At the two Reynolds numbers the starting flow is similar except that 
vorticity is spread over a larger area for R = 40 than for R= 100. ‘The 
velocity within and near the rear eddies is also much smaller for the lower 
Reynolds number. For the drag coefficients in the two cases see figure 1. 
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The mean velocity of slightly buoyant and heavy 
particles in turbulent flow in a pipe 
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SUMMARY 

A large number of small spheres of the same size were injected 
successively into a horizontal pipe conveying water at constant 
mean velocity, and their times of transit were measured. ‘The 
mean velocity of the spheres that were either somewhat heavier or 
lighter than water was less than that of those of neutral density ; 
for those having a terminal velocity in water within +1°,, of the 
mean velocity of the water in the pipe, the discrepancy was only 
about 0-1°,,.. The dispersion of the times of transit of the spheres 
was almost independent of their density. A theory is developed 
to show how the mean velocity of the spheres depends upon their 
relative density and size. 


1. INTRODUCTION 

In a previous paper (Batchelor, Binnie & Phillips 1955) an account 
was given of experiments in which numerous spheres of almost neutral 
density were injected in turn into a turbulent stream of water in a horizontal 
pipe and their passage timed. In this way a method was devised for finding 
the discharge of water through the pipe and the fluctuation in the times 
of travel of the spheres. It was shown theoretically that the mean velocity 
of material particles of the liquid is equal to the discharge velocity U, 
defined as the discharge, at some cross-section, averaged over a long time 
and divided by the cross-section. A satisfactory modification was developed 
to account for the difference between U and the mean velocity U(x) of spheres 
of small but finite size, this difference being a function of x, the ratio of the 
diameter of the spheres to that of the pipe. ‘The density of the spheres 
was such that their terminal velocity in water was within +1°,, of the 
discharge velocity, as determined by a measuring tank at the pipe outlet. 
‘This standard of density is not easy to achieve under the conditions of 
temperature and pressure prevailing in the pipe; and before the method 
can be used with confidence to determine the discharge velocity through 
a large pipe-line, it is important to know the consequences when the above 
standard is slightly relaxed. An attempt to answer this question by means 

At present at Mechanical Engineering Department, The Johns Hopkins 
University. 
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of further experiments at the Engineering Laboratory, Cambridge, is 
described below. 

If the spheres are buoyant or heavy, their motion is modified in two 
ways. Firstly, the difference between the inertia of a sphere and that of 
an equal volume of fluid results in relative motion when the fluid surrounding 
the sphere is accelerated. Secondly, the effect of gravity is to modify the 
probability density of the sphere positions in such a way that for spheres 
denser than water the probability density is greater towards the bottom 
of the horizontal pipe than towards the top, and vice versa for spheres 
lighter than water. It may reasonably be supposed that the effect of gravity 
on the motion of the spheres is determined by the value of y, the ratio of 
their terminal velocity V to the discharge velocity U. The results of the 
experiments described in §3 suggest strongly that the inertia effect is 
negligible when y is small but that the gravity effect is not; we shall 
anticipate this result by writing the mean velocity of the spheres as U(z, y). 

In the light of experience the apparatus previously used was rebuilt 
in the manner explained in §2. ‘The most important improvements were 
the introduction of an automatic device for injecting the spheres and of 
better methods of controlling and timing the discharge. More robust 
spheres were required in the new injection apparatus, and they were made 
of polythene in place of wax. In the experiments attention was confined 
to spheres of diameter 0-2 in. in a pipe of diameter 2 in., hence « was fixed 
at 0-1. ‘The Reynolds number, based on pipe diameter, was about 7 x 10+. 

The experimental results make it clear that, to achieve the highest 
accuracy, the 1°,, standard of terminal velocity should be adhered to. 
However, the available evidence suggests that little advantage will follow 
from efforts to improve on this standard. Finally, in §4, an expression for 
U(z, y) is derived theoretically, and with the aid of the experimental results 
tentative predictions are made of its magnitude for values of « other 
than 0-1. 


2. DESCRIPTION OF APPARATUS 

‘The entire pipe-line was rebuilt in new galvanized steel tubing, the mean 
internal diameter of which in the working section was found by weighing 
the water required to fill each length. ‘The supply was through a 6-in. main 
from a large tank about 40 ft. above, in which the level was kept constant 
within +1in. Near the inlet a stainless-steel orifice-plate was inserted 
to act as a resistance additional to the wall friction. ‘The experiments were 
nominally conducted at one velocity only, and all that was required to 
control the water flow was the complete opening of a valve. Nearby, the 
spheres were injected by means of the device shown in elevation in figure 1. 
The long tube A, which extended to the centre-line of the pipe, could 
carry a charge of 100 spheres. ‘The spheres in it were pushed downwards 
by a small stream of water, supplied at a steady net head of a few inches 
to the top of the tube and escaping at the bottom into the pipe. The electric 


motor B revolved the disc C, which near its periphery carried the vertical 
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pin D. Close and parallel to the tube A a shaft was mounted to which a 
small arm was fixed at each end. ‘The upper arm protruded above the 
disc, so that every 15 sec the shaft was rotated by the pin D through 90 
and was then returned to its original position by a spring. The bottom 
arm, which is shown more clearly in the plan and elevation drawn to a 
larger scale on the left of figure 1, engaged in a slot the lowest sphere in 
tube A; thus, when the shaft rotated, the sphere was moved out into the 
main stream. Release was assisted by a small vane that deflected the stream 


downwards into the slot. 
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Figure 1. Injection apparatus. 


After passing round a bend the water entered on a long, straight and 
horizontal course. The first 28 ft. served to remove the disturbance due 
to the bend. The remainder was used for the timing measurements, and 
in it three lengths of Perspex tube were inserted, fitted with the same 
arrangement of three photo-cells and two Dekatron counters as was 
previously employed. Each sphere at the first station A started both 
counters, and stopped them in turn as it passed stations B and C. The 
lengths AB and BC were respectively 33-42 and 16-40 ft., and each 
injection provided two observations directly and one by difference. At 
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exit the water passed into a pool where the spheres were netted and then 
over a weir into a measuring tank. ‘There the rise was timed at regular 
intervals with the aid of two vertical probes which started and stopped 
another Dekatron counter. All the counters were operated by the mains, 
and long-period errors in the mains frequency did not atfect the measured 
ratio of the velocity of the sphere to the discharge velocity. 

Spheres made of the wax previously used were found to be too soft 
when tried in the new injection device, and this material was replaced by 
sheet polythene loaded to have a relative density close to unity. At first, 
the spheres were made in a cold press and by hand moulding after the 
material had been softened by heat. ‘These processes accidentally produced 
considerable numbers of spheres in which a small air bubble was trapped. 
Later, hot injection into a mould was employed usually yielding slightly 
heavier spheres, and some of these were made heavier still by driving in 
the point of a pin. ‘The terminal velocity of each sphere was measured 
in a tube 58 in. long. and the spheres were batched in terms of y, the ratio 
of this velocity to the discharge velocity, which was about 5-37 ft./sec. 
The head in the working section of the pipe was about 3 ft., and the spheres 
were tested at this pressure in case the presence of the bubbles made the 
spheres appreciably compressible, but in fact the terminal velocity even 
of the lightest spheres was found to be independent of the depth of immersion 
in the tube. ‘The temperature coefficient of expansion of polythene is 
greater than that of water, and the temperature of the water in the laboratory 
circulating system was not under control; but its variations were slow, 
and each batch of spheres was checked at the correct temperature immediately 
before use. Sufficient spheres were made to permit observations with five 
batches, for which y lay within the limits 0-04 to 0-03, 0-02 to 0-01, 0-01 to 0, 

0-01 to —0-02, —0-02 to — 0-03, and a further batch at —0-0375, which 
was the mean for a mixture of batches — 0-03 to — 0-04 and — 0-04 to — 0-05. 
Here the negative sign denotes spheres lighter than water. At 15°C the 
relative densities of spheres for which y = 0-04 and 0-01 were calculated 
to be 1-042 and 1-005 respectively. 


3. EXPERIMENTAL RESULTS 


f y, 200 observations were made of the times of transit 


For each value of ; 
of the spheres between the cross-sections 4B and AC, measurements of 
the discharge being taken at the same time at regular intervals. “he mean 
velocity U(z,y) of the spheres was determined from equation (4) of the 
previous paper, which is 

U(a,y) = x/T(x), 
where 7(x) is the average time taken for the spheres to travel a length x 
of the pipe. In this equation there is a proportional error due to dispersion, 
which can be calculated with the aid of equation (19) of the previous paper. 
For x = 0-1 it amounted to about 2a/x, where a is the pipe radius; and it 


was ignored because even for the shortest length BC it was only 0-1°,. 
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The measured mean diameters of the lengths AC, AB and BC were 
respectively 2-005, 2-009 and 1-998 in. (the variation of diameter along each 
length being unknown). ‘The first of these values was employed in 
calculating U from the discharge measurements, and the results are shown 
in figure 2 in the form {U(a,y)— U}/U vs y. Now figure 2 of the previous 
paper shows that for spheres of the same density as water dj U(x)/U}/dx 
is roughly } near « = 0-1. So, if the pipe cross-section is diminished by 
1°, thus raising U by the same amount and « by }°%%, the increase in 
|U(a)— U}/U is 0-0025. The correction to the observations in AB and AC 
is therefore small; and it has not been made in figure 2 since a greater 
precision in the absolute value of U cannot be claimed. 
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Figure 2. Effect of density upon mean velocity. eAC, x AB, + BC. 


Near y = 0 the results lie somewhat below the value 0-075 that was 
found in the earlier work for « = 0-1. The discrepancy is probably due to 
the relatively crude method then used for determining the discharge. 
Even with the more refined method used in the present work, it is difficult 
to determine U with the accuracy that 1s readily achieved in the measurement 
of U(a,y). 

‘The variation of the mean velocity of the spheres with relative density 
may be due to inertia or to gravity effects (or to a combination of both). 
The former effect is difficult to analyse, but it is unlikely to be the same 
for equal and opposite departures of the relative density from 1:0. The 
gravity effect, on the other hand, is obviously the same for light and heavy 
spheres whose densities differ from that of water by amounts of the same 
magnitude. It is also probable that the effect of gravity is to diminish 
the mean velocity of non-neutral spheres, because the trajectory of lighter 
or heavier spheres is to a greater extent in regions near the wall at the top 








92 A. M. Binnie and O. M. Phillips 


or bottom respectively of the pipe, where the velocity of the water is 
smaller. We do in fact see from figure 2 that the results are almost 
symmetrical about the line y = 0, therefore it was probably the effect 
of gravity which produced the changes in U(z,y) as y was varied. 
The experimental curve is fairly flat over the range 0-01 > y > —0-01, 
but falls rapidly outside these limits. Thus, for « = 0-1 at any rate. 
spheres made within these limits will move with very nearly the same 
mean velocity as spheres of precisely neutral density, the error being no 
more than 0:1°,,. 

For a pipe inclined at an angle ¢ to the horizontal, figure 2 might perhaps 
be used with y taken as (Vs U)cosd. However, if ¢ is not a small angle 
it may be necessary to allow for the gravitational drift of the spheres along 
the pipe due to the component V sing. 
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Figure 3. Effect of density upon the dispersion coefficient. ed4AC, x AB, + BC 


Figure 3 shows the values of the longitudinal dispersion coefficient 
K(z,y), defined as in equation (19) of the previous paper by 
v U [T(x)- T(x)}? 


K(x, y) 


/ 


2a u. T( v) 
where w, is the friction velocity, and 7(x) is the time of travel of a single 
sphere. Unlike {U(z,y)—U}/U, the value of A is not sensitive to small 
changes in y, and the agreement at y = 0) with the earlier work is satisfactory. 
3ut the accuracy of these measurements involving the variance of the travel 
times is not high. ‘The standard deviation of the measured value of K is 
given by elementary statistical theory as a fraction n~!? of the true value, 


where 7 is the number of observations. Here 2 = 200, so that the expected 
standard deviation is about 14°,, which is in accordance with the scatter 
found. However, the results do seem to have a trace of a minimum near 


y=. Since the longitudinal dispersion is largely a result of the 











Mean velocity of particles in turbulent flow 93 


intermittent excursions of the spheres into the region of lower velocity near 
the wall, an increase in K(x, y) when y is different from zero is consistent 
with the above interpretation of the results shown in figure 2. 


4, ‘THEORETICAL DISCUSSION 

There remains the question of accounting for these observations by an 
approximate theory so that predictions can be made concerning the 
behaviour of {U(«,y)—U}/U for sphere sizes other than that (« = 0-1) 
used in these experiments. If the spheres are not of neutral density, that 
is, if y # 0, the probability density p of the sphere distribution (defined so 
that the probability of finding the centre of the sphere in an element dA 
of the cross-section is pdA) is not uniform over the cross-section. ‘The 
observations show that inertia effects can be neglected as a first approximation 
when y is small, and it can be supposed that as a result of gravity the 
probability density p is a function of the vertical position coordinate 2 
only, that is, that the lines of equal probability density are horizontal. 
Then, if V represents the terminal velocity of the spheres (a positive sign 
indicating motion upwards) and if conditions are statistically stationary 
in time, the probable number of spheres rising across unit horizontal area 
per unit time as a result of the action of gravity alone is p(z)V’, which must 
be balanced by the downward flux due to turbulent diffusion, u, a dp(z)/dz, 
where uw, is the friction velocity and ¢ is a dimensionless turbulent diffusivity 
for transport across horizontal planes. In general, ¢ is also a function of 
position in the pipe, and depends upon the relative intensities of the 
turbulent fluctuations in the vertical direction at different points in the 
cross-section. However, if the sphere size is not so small that its motion 
is affected by excursions into the viscous sub-layer, it is confined to the 
central portion of the flow where the turbulent intensities do not vary by 
a factor of more than about 1-5. A further approximation, consistent with 
taking p = p(s), is to neglect the variation in ¢ with position, and to consider 
it as a mean vertical diffusivity for the turbulent flow in a pipe. We then 








have ' 
dp(z I 
a: id SO (1) 
dz u_ac 
so that ; pew 
plz) = poexp(yl s"/u,C), (2) 
where z’=s/a. The constant of integration py is determined by the 


condition that p(z’) integrated over that part of the cross-section accessible 

to the centre of the sphere, that is, a circle of radius (1—«)a, must equal 

unity. When y is small, the exponential can be expanded as a power series, 
and the integration is found to give the result 

¢ a i ¥\9 2 

po = (7a?(1 —)} 1 — 3 (yU/u, CP? (1-2)? + O(y*)}. (3) 

Let us now make the further assumption that the mean velocity of 

spheres whose centres are at a radial position r in the pipe is equal to the 

mean velocity u(r) of the fluid at radius r, which was shown in figure 2 of 
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the previous paper to be a good approximation for « < 0-15. ‘Taking 
zx = rcos# andr’ = r/a = z'/cos§, we see that the mean velocity of a sphere 
down the pipe is given by 


‘2x pl—o@ 


U(a,y)=a?] | p(s’ u(r’ )r’ dr'dé. (+) 


Substituting from (2) and (3) into (4) and integrating with respect to 4, 


we have, neglecting terms of order y', 
. . vl/ 2 ‘l~« pia 
U(z,y) = U(x)- (5) <i(1—a)? | ru(r)dr—$ | = ru(r)dr>, 
(5) 


where the accents on r’ are suppressed and 
“22 -l-—a@ 
U(a) = U(a,y = 0) = {x(1—«)*}-! | | u(r)r drdé. 


Taylor (1954) shows that u(r) is of the form 

u(r) = U+u{4-25—-f(r)}, (6) 

where f(r) isa universal function, and that for Reynolds numbers of order 10”, 
as in our experiments, U/u, = 22. Hence 





u(r) = U{1-193 —0-0455 f(r)}. (7) 
Equation (5) can then be expressed as 
SMe) 7 ON) _ _ Fl -@)'Ky@)-2K@), 8) 
( (1—a)? & 
where 
-l—« -1—x 7 
K, (a) ji | ru(r) dr = 0-597(1 —«)* —0-0455 | rf(r) dr, 
l al ns y (9 
K,(«) = U | ru(r) dr = 0-299(1 —«)*—0-0455 | rf(r) dr. 
Jo Jo J 





x | K,(a) K;(a) 


0-02 | 0:-4899 | 0-2179 
0-04 | 0-4749 00-2046 
0-06 0-4599 0-1913 
0-08 0-4460 0O-1781 
0-10 | 0-4341 | 0-1650 








lable 1. Values of Ky(x) and K,(a). 


The integrals on the right of the equations (9) were calculated from 
the figures given in ‘T'aylor’s table 1, and the quantities A,(x) and A3(z) 
are given in table 1 for various small values of x. For « = 0-1 it is 
then found from (8) that 


U(a,y)— U(a) m ey a (10) 
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‘The experiments with x = 0-1 described in §3 enable an estimate to be 
made of the dimensionless diffusivity ¢. The ratio {U(x%,y)—U}/U for 
z = 0-1 and y = 0 was judged to be 0-0615 from figure 2; and by subtraction, 
values of | U(x, y)—U(x)}/U were obtained. The parameter ¢ was then 
chosen so that the parabolic form (10) fitted these observations most closely, 
and the value required was found to be 0-46. ‘The mean vertical diffusivity 
is therefore 0-46au,. It is interesting to notice that the radial eddy 
diffusivity in turbulent pipe flow found in another way is of the same order 
as this, some experiments of Schwarz & Hoelscher (1956) with water vapour 
in air indicating a mean value of approximately 0-55au, in the central 
region of the pipe. 
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Figure 4+. ‘Theoretical effect of density upon mean velocity. 


With the value of € determined above, the form of {U(z,y)— U(«)}/U 
for small values of y and for x = 0-06 and 0-02 was found from equation (8), 
and the results are shown in figure 4. For a given value of y, the difference 
between the mean velocity of the spheres down the pipe and the mean 
velocity of those for which y = 0 increases with decreasing diameter ratio x. 
The physical reason for this greater etfect of gravity on smaller spheres 
is that small spheres are free to move closer to the walls where the mean 
velocity of the fluid is less and is changing rapidly; the net effect, then, of 
the increased probability of finding (say) a heavy sphere near the lower 
wall and the decreased probability of finding it near the wall in the upper 
part of the pipe is that, for a given value of y, (U(x, y)— U(«)}/U increases 
in absolute value as ~ decreases. ‘hese curves indicate that, provided y is 
held within the limits +0-01, the difference between the mean velocity 
of such spheres and of those for which y is accurately zero is likely to 
be negligible, as has already been established experimentally for « = 0-1. 
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We are greatly indebted to the Plastics Division of Imperial Chemical 
Industries Ltd., who gave us the loaded polythene and carried out the 
hot moulding of the spheres. 
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Water waves of finite amplitude on a sloping beach 
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SUMMARY 

In this paper, we investigate the behaviour of a wave as it 
climbs a sloping beach. Explicit solutions of the equations of 
the non-linear inviscid shallow-water theory are obtained for 
several physically interesting wave-forms. In particular it is 
shown that waves can climb a sloping beach without breaking. 
Formulae for the motions of the instantaneous shoreline as well 
as the time histories of specific wave-forms are presented. 


1. INTRODUCTION 

The behaviour of waves on sloping beaches has received intensive study 
by many authors during the past sixty years. ‘These investigations, for the 
most part, have been confined to studies of linearized problems which are 
based on assumptions that are invalid in the neighbourhood of the coastline. 
With the results of these linear theories as a basis, it has been stated that 
progressing waves eventually break on a sloping beach. 

In this paper, we present an analysis based on the non-linear shallow- 
water theory. Explicit solutions are obtained for a number of important 
cases and, in particular, it is shown that there are waves that climb a sloping 
beach without breaking. ‘The initial shape and particle velocity distribution 
determine whether or not a given wave will break, and no simple criterion 
for the occurrence of breaking has been found. 


2. GENERAL ANALYSIS 
The conservation equations of mass and momentum of the non-linear 
shallow-water theory are 


[o*(n* +h*)] ae = — Tye (2.1) 

vs t+vto%, = —Z7%, (2.2) 

where 7* is the horizontal velocity and the other symbols are defined in 
figure 1; the asterisks denote dimensional quantities. A complete 


development of the non-linear shallow-water theory can be found in Stoker 
(1948). 

It is convenient to introduce the following dimensionless quantities : 
v= vl, n= 7*/ah, x= x*/l, t=t/T, & = (h*+7P)/al. In these 
definitions, 7 = (j/xg) and vy =(glx)'?. ‘The characteristic length /, 
can be specified when a specific problem is adopted for study; the depth 
is assumed to be of uniform slope, * = —ax*, 
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With these substitutions, equations (2.1) and (2.2) become 
v.t+e¢v,+n, = 0, (2.3) 
) 


[7(m—x)],.+ 7, = 0. (2.4 
‘These hyperbolic equations can be rewritten in a form in which the 
characteristic variables x, 8 play the role of the independent variables and 
7, ¢, x and ¢ play the role of the unknown functions of «, 8. ‘The four 
equations which arise when the classical transformation is made (the 
details are given in Stoker (1948)) are 


v,—(v+c)t, = 0, x, —(v—c)t, = 0, (2.5), (2.6) 
Ugt+2eg+t, = 0, U,—2¢, +t, = 0. (2.7), (2.8) 

Equations (2.7) and (2.8) can be integrated explicitly to obtain 
v+2c+t=4, vw—2c+t= —B. (2.9), (2.10) 








Pg 
Figure 1. Definition sketch. ‘The fluid has a sloping fixed boundary and a free surface 
at height 7* above its undisturbed level. 
Here, the ‘constants of integration’ have been chosen in the interest of 
algebraic simplicity in what follows. From (2.9) and (2.10) we obtain 
v+t =(a—B)/2 = A/2, (2.11) 
c = (a+ B)/4 = o/4; (2.12) 


and these define A and o. We now adopt 4 and o as our final pair of 
independent variables, so that (2.5) and (2.6) become 


v,—wvt,+ct, = 0, (2.13) 
v,+ct,—vt, = 0. (2.14) 

The elimination of v results in the /imear second-order equation for f 
ai,,—t,.)—3t, = 0; (2.15) 


and, since 7 +1 = A/2, v must also satisfy (2.15). In fact, it is readily verified 
by reference to equations (2.11) to (2.14) that if we introduce the ‘ potential’ 
d(o. 2) so that 


uo a '¢,(a, A), (2.16) 
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then 
x = $,/4-—07/16 —v2/2, (2.17) 
n= C+x = 6,/4—-v?/2, (2.18) 
t = A/2-2, (2.19) 

and 
(o4,,), —06;, = 0. (2-2) 


Two major simplifications have been obtained. ‘lhe non-linear set 
of equations (2.5) through (2.8) have been reduced to a linear equation 
for v or ¢ and the free boundary (the instantaneous shoreline c = 0, which 
moves as a wave climbs a beach) is now the fixed line o = 0 in the (o, A)-plane. 

The choice of a function ¢(¢,A) which satisfies equation (2.20) defines 
n, U, xX, t in terms of the parametric coordinates o, A. In particular, if the 
Jacobian 0(x,t)/e(o,A) never vanishes in o > 0, the implicitly defined 
solutions 7(x,¢) and v(x, t) are single-valued, and such solutions represent 
waves which do not break. If the foregoing Jacobian does vanish in o > 0, 
the wave must break. However, we confine our attention in this paper 
to those forms of ¢ for which breaking does not occur. 

A particularly simple solution of these equations is given by 


6 = AJ ,(wa)cos(wA— ys), (2.21) 


where J, is the usual notation for a Bessel function. No loss in generality 
ensues when the phase lag 7 is taken to be zero or when w& is put equal 
to unity, so that 

6 = AJ,(a)cosa. 


The Jacobian J = o(x,t)/e(o,A) vanishes nowhere in o > 0 when 
A <1 and the mapping is valid in o > 0. ‘The physical phenomenon 
whose description is implied by equation (2.21) is that which occurs when 
a wave of unit frequency in the dimensionless time variable travels 
shoreward from the region of large x and is reflected, so that a wave, again 
of unit frequency, travels out to sea. ‘The reflection coefficient is unity. 
The phenomenon is periodic in the time variable and the wave shape 
far at sea is like J,(4\/ x ), but it is considerably distorted near the shore. 
In particular, the penetration of the wave (the value of w at which the depth 
is zero) is given by equation (2.17) when o = 0: Le. x(A,0) = 6,/4—u?/2, 
so the maximum penetration is 4/4. When 4 > 1, J vanishes on some 
curve ino > 0, and the solution must be modified so that a bore is included 
in the prediction. ‘The analysis of that problem is now being studied, 
but will not be discussed here. ‘The wave shape is shown in figures 2 and 3, 
for the extreme positions corresponding to A = 7/2 and A = 32/2 both for 
A=1and A=}. In the limit as 4 —-- 0, J9(a) becomes J,(4y/ x ), the 
linearized solution, and no graph should be needed. 

We now consider problems in which a mound of water is released: 
that is, we specify a wave shape 7(x, 0) with v(x, 0) = 0 everywhere. Since 
7+¢= 4/2, the condition that 7 = 0 when t = 0 implies that A = 0 for 
t=. Equation (2.15) for 7 must be solved in the region o > 0, A > 0, 
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with the initial conditions that () and @, is specified on A = 0, and the 


boundary condition that v is to be finite on a = 0. 
The derivative v, can be determined from the prescribed initial wave 

height by first solving the equation 
[4(x, 0) —x]#? = e(x, 0) = 


for x as a function of o, and then using equation (2.13) to show that on A = 0 


] 
1o 


x. ct.. 


Figure 2. Free surface geometry for periodic motion according to equation (2:21) 
for A |! at: E, point of maximum penetration, A 7/2; B, A Saiz: 
(’, D, intermediate times. 
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| 
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Figure 3. Free surface for periodic motion according to equation (2.21) for A 
E, a 7/2; B, aA 37/2: ©, D, intermediate times. 


at: 


‘Therefore, for A (), 
7 I ip + 4oq 2 f(o). 


\n entirely equivalent boundary-value problem can be formulated for the 


potential ¢. 
The solution of the boundary-value problem for the particle velocity 7 


is easily obtained by transform techniques. Let 
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then equation (2.15) becomes 


ov 4 3 U_—s“ov = of(o). (2 Ze 
Let z = so and p = 20; then 
(xp’) ma 2p = — — f(2/s), (2.23) 
g s* 


and if we now use a Hankel transform 


rw 


p= | oJ (Ez)p(z) dz, 


we obtain after further manipulations 


. 1 (4 
P- TrB. 


ke 
t 


J (Bs) f(B/s) dB. 





9 


gy § 


Upon setting € = st and computing the inverse transforms, there results 


U = | a (ra)sin tA dr | a2 J ,(ray) f(a) doy, (2.24) 
or, in terms of the potential, 


o=— | tJ (7ra)sin tA dr [ o- J (704) f(a) doy. (2:25) 


ae |) 0 


These results can also be obtained by the superposition of standing 
waves. The function o1J,(ko)sinkA is a solution of (2.21). By the 
principle of superposition, 

v= | A(k)oJ,(ko)sin kA dk 
is also a solution. The function A(k) is determined from the boundary 
condition v, = f(c) onA = 0; the condition that v = 0 on A = 0 is implicitly 
satisfied. Further reduction of these general expressions leaving f(c) 
unspecified does not simplify the task of evaluating the final integrals. 
Instead, we select functions f(a) which will both simplify the final integrals 
and correspond to physically interesting initial wave shapes. 


3. SOME INITIAL VALUE PROBLEMS 
In this section we consider a number of interesting examples of wave 
propagation problems in which the motion starts from rest at time zero. 
As a first example, let a one-parameter family of wave-forms at t = 0) be 


given by 





o 5 a 3 a 
v= —— +e 1- >752 232 «0/2 2 
16 2(a*+0*)®? 2(a*%+0°*)* 


where a = 1:5(1+0-9e)'!*. These waves, shown in figure 4+ for two values 
of ¢, all have maximum heights equal to ¢, all have heights 0-9e at x : 
and all have zero slopes at the shoreline. ‘These shapes correspond to the 
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fluid held motionless and then released. 


equ 


and 


Figure 4. 
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sical problem in which the water level at the coastline is depressed, the 
‘The quantity f(c) is found to 


} 


| " a2 J ,(70y)f (a9) doy = 2a®ere~*" (4 — ar). 


al 
1 


+ g?)>2 


2 
30a%e ; : — 
(a (a? + 0”)? 


for this function it can be shown that 
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(3. 
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Initial wave shapes given by equations (3.1) and (3.2) for «-- UO and 


0-1. 


€= 


If we solve for v using equation (2.19), set o = ao’, A = aX’ and then 


drop the prime notation, we find that 


and 
‘The 
obti 








8 ; 1 3 1—2A 
oo ie a — ———- |. (3.4) 
a (1 —2A)? +0738" = -4:4(1 —21A)? +07)? 
1 1 1—1A ‘ 
sed Ion eee: —F— ae - = eT (3.5) 
'(1—2A)®?+o2}2  4/(1—i1A)? +0?}8? 
Uv a’o? a > 1 2 3/4 1A 3 ate! ll = 3 (3.6) 
2 16 f(1 —2A)?+07}82 2 {(1 —2A)? + 0?}? 
= la\—v, c = jaa, (3.7), (3.8) 
y = x +a?o?/16. (3.9) 
* motion of the instantaneous shoreline, or zero depth position, 1s 
ned by setting o = 0, and is given by 
>. £3 _ 5 
Se SA" ~ (3.10) 
a (1+A?) 
‘ 
: qe 4 (ae), eee oye 
x € (i | ) ( ) 
t= la\-i c= 0. (3.42), (3.13) 
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The maximum penetration distance attained by the climbing wave 
occurs when the coastline velocity is zero, i.e. for A* = 5, and is given by 
Xyay = 1:157e. ‘That is, the water level, if depressed a depth ex/) at the 
shoreline and then released, will rise to a height 15°, higher than the 
original sea level. Figures 5, 6 and 7 present a time history of the action, 
and figure 8 is a plot of the position and velocity of the instantaneous 
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Figure 5. Time history of the wave-form of equation (3.1) for « = 0-2, near the 
coastline. 


shoreline for the specitic wave shape e = 0-1. Itis seen that the instantaneous 
shoreline rises above the mean sea level and then slowly settles back to it. 
There are no continued oscillations about this position and the waves 
do not break provided ¢ is sufficiently small, namely « < 0-23. 
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As a second example consider the motion of a stationary mound of 
water released at ¢ = 0, which is given by 
n = tepreate—?, (3.14) 
and 
x = Jee*p*ote °? — 07/16, (3.15) 
where 1/p = 8(1+.e). ‘The family of wave-forms, shown in figure 9, all 
have their initial maximum wave heights at a fixed position from the shoreline, 
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Figure 6. Time history of the wave-form of equation (3.1) for ():2, far from the 
coastline. 


7» = eatx = —1. Inaddition all have zero slope at the origin. The quantity 
v, evaluated at A = 0 is found to be 
f(a) = 2ee*p?(207e-7*” — ope”) ; (3.16) 


and from Watson (1944, § 13.3) we find that 


FP a es e. nf | eo 
| or * J (ray )e P doy (-1) - (pe *) (3.17) 
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MAXIMUM WAVE SLOPE 
7), ~100 








Figure 7. ‘lime history of the wave-form of equation (3.1), € 0-2, in the neighbour- 


hood of the beach. 











! \ 
eb fee ee 
r eae —____ 
yr | if 
4 \ 
om /\ 1.08 
| / 
= i \ 
| \ = 
4 | | \ 41.04 
Bie a 
| | ia \ 
a / \ — 1.02 
4 \ 
r/ ’ a 
/ it 
4 8 2 € OO er Oe 5.0 t 
S na 











Figure 8. Coastline position and velocity vs time for the wave-form of equation (3.1) 
with « = 0:1. 
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If we set s = itp '*, then the complete solution is 
r tee- | SO ' (2p! *az)sin(2p! *As)e "Ee 2: t ! 24) dz, (3.18) 
hb = Zee*p * | sin(2p'Az)J,(2p!a2)e *(1 — 22? + Is) dz, (3.19) 
v -v?/2 —o7/16+€e* | 3 cos(2p! *A)e-* (1 — 23? + }34)J,(2p! 202) dz, 
(3.20) 
t= A/2—9¢,/¢, and c= fo. (3.21), (3.22) 
We 
i i 'S ye ie 9 
L € i 
i 1 i 1 1 4. 
€2#j5 4 
R e=| ad 
4, 4 1 4 
8 7 ~€ rs -4 -3 2 - = 6) 
Figure 9. Exponential wave-forms of equation (3.15); « 0, 0°10, €. 


[t is seen that the quantities 7, 7 
independent of o and A: Le. 
v| < Me, v,| < Mie, v,| < Mge. 


lor € sufficiently small, the Jacobian 


and v, have upper bounds which are 


4 


J|=1e v2—(! —wv,)* 
is approximately given by o/16, and therefore does not vanish in the interior 
of the fluid. In other words, for « small enough, the wave forms of 
equations (3.14) and (3.15) do not break as they climb the beach. That 
the Jacobian is zero when a is zero is a property of the transformations and 
is not a consequence of assuming a particular initial wave shape. The 
slope of the wave at the instantaneous shoreline, o = 0, can remain finite 
even though the Jacobian vanishes. Similar statements hold .or the wave 


shapes of the previous example. 
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If we restrict ourselves to a discussion of the motion of the shoreline, 
set A’ = 2p'*A, o’ = 2p'*o and afterwards drop the prime notation, then 
the equations of motion of the shoreline can be written as 


U (zp)! =e? df (A) da, (3.23) 
t = }Ap-"? — (ap)! ee? df(A)/da (3.24) 
c=0, x = —v*/16 + €e?7!? f(A)/4, (3.25), (3.26) 


where 











a 
Figure 10. A plot of the functions 16f and 16f, ws A. 


‘The functions 16f, 16f, are shown in figure 10; asymptotically 16f ~ — 32/A*, 
lof, ~ 64/A, so that v ~ 4(zp)!ee?/A®. "The coastline motion is such that 
it first rises to its maximum height, falls back below its initial position, and 
finally returns, very slowly, to the original mean sea level. ‘There are no 
continued oscillations about the mean sea level. ‘The maximum penetration 
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distance occurs when wv = 0, for A = 2-41 and is found to be x,,,, = 1-45le. 
‘The maximum height attained is 45°,, greater than the maximum initial 
wave height. The coastline velocity is again zero at A = 4-835, which implies 
that the lowest depth reached is x,,,, = —0-636e. ‘The shoreline motion 


for the particular wave shapes € = 0-1 and ¢€ = 0-5 are shown in figures 11 


and 12. 














; 


Figure 11. Coastline position and velocity vs time for the exponential wave 
e = 0-1 of equation (3.15). 


4. CONCLUSION 
‘Thus far we have shown that there are progressing waves of a com- 
pressive nature (positive amplitude) which do not break as they climb a 
sloping beach. No general criteria have been found which enable us to 
determine when a given wave will break, although the magnitude of « 
and hence the initial wave shape are of fundamental importance. The 
wave-forms we have considered all had continuous derivatives, the first 


derivative being zero, initially, at the coastline. Such waves do not break 
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| 
Figure 12. Coastline position and velocity vs time for the exponential wave 
€ 0-5 of equation (3.15). 


for sufficiently small «. In a subsequent paper, it will be shown that all 
compressive waves (waves of positive amplitude) propagating into 
quiescent water towards the beach with a discontinuity in first derivative 
necessarily break before reaching the coastline. 
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REVIEWS 


An Introduction to Fluid Mechanics and Heat Transfer, by J. M. Kay. 
Cambridge University Press, 1957. 309 pp. 37s. 6d. or $7.00. 


The development of a satisfactory syllabus for students of engineering 
at the university is a difheult undertaking. Particularly is this true in the 
teaching of chemical engineers, the jacks-of-all-trades, or where the intention 
is to give a general engineering training without specialization. Professor 
Kay’s book, which owes much to his experience of instructing chemical 
engineers at Cambridge University, clearly reveals the difficulties of compiling 
a course in one of the subjects that engineers study and the ingenious efforts 
which have gone towards overcoming these difficulties. 

‘The problems that beset the university engineering instructor arise 
from a dilemma: must he make a choice between presenting an academic 
discipline and giving vocational training, or is a workable compromise 
possible? ‘There is too much pressure from other subjects in the curriculum 
to allow full formal development of topics, nor would the mathematical 
equipment of the average engineering undergraduate permit it. Quite 
understandably, most student engineers are too eager to be applying their 
knowledge to be ready to spare much time and attention for the finer points. 
Kay definitely chooses to attempt the compromise. His book aims to be 
fundamental and systematic enough to provide some intellectual satisfaction, 
but the utility and possible applications of the analysis are emphasized at 
every stage in order to match the engineer’s constructive outlook. 

More than most, fluid dynamics is a difficult subject to teach to under- 
graduate engineers because only a limited amount of material can be taught 
with the aid of elementary mathematics of the standard also needed for the 
courses on mechanics, structures, thermodynamics and electricity, that the 
student is usually taking at the same time. ‘This mathematical limit 
restricts the fluid dynamics course to what can be done by means of simple 
momentum statements, Bernoulli’s equation and dimensional analysis, 
bolstered with empirical and theoretical facts taken on trust. ‘That under- 
graduates’ fluid dynamical training should often stop at this elementary 
level is very regrettable in view of the dominance of problems involving 
Huid motion in many current scientific and technological activities, as the 
pages of this journal show. Indeed, most newly graduated engineers would 
tind the reading of this and kindred journals well beyond their immediate 
capacity. 

Kay is evidently eager to remedy this situation, but naturally the trans- 
formation of undergraduate curricula cannot be achieved in one step. ‘The 
result is what one might call a transitional book, where the author breaks 
from the customary limitations, but does not entirely find a satisfactory 
alternative. Vector notation must be adopted, the author decrees, and one 
applauds the decision, but still much of the essential matter is stowed away, 


half apologetically, in appendices. And if vectors are to be used, why should 
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not elementary cartesian tensor notation also be introduced, particularly 
in the discussions of stress, strain and vorticity ? 

The separation of fluid mechanics from thermodynamics has become 
increasingly artificial and one may expect the appearance of more and more 
books which straddle the two subjects as this one does. As its title would 
suggest, the book contains only thermodynamical topics of relevance to 
fluid mechanics and heat transfer. ‘The thermodynamic fundamentals are 
not well developed. One could wish for a stricter treatment of thermo- 
dynamical concepts and the abandonment of the term frictional heat where in 
fact no heat is flowing. ‘he Second Law is neglected to a surprising degree. 
It is nowhere explicitly stated, although the significance of the sign of the 
entropy change in shocks is mentioned briefly. It is pleasing to find the 
parallelism between momentum, heat and mass transfer explained well, but 
there is need for a more forceful statement of the one-way nature of the 
interaction between the fluid motion and the heat transport in low-speed 
forced convection if the effect of temperature on fluid properties is negligible. 

In view of Professor Kay’s interest in the field of nuclear power, where 
fluid dynamicsand heat transfer are closely inter-related, itis surprising to find 
so few references to nuclear engineering problems. _‘I’o choose an instructive 
example, the inevitable concomitance of wall friction with heat transfer in 
pipe-flow is the dominant factor in the problem of cooling power reactors 
without excessive degradation of hard-won power in circulating the coolant. 
Heat transfer in high speed gas flow is another topic, important in many 
fields, which Kay’s book neglects, while heat transfer at very low Prandtl 
numbers with liquid metals also receives scant attention. 

The book is evidently based on a lecture course which occupied more than 
one year. One symptom of this is the way that topics such as the normal shock 
tend to reappear in more sophisticated detail as the book proceeds. Another 
is that facts are often asserted in the earlier chapters with references to proofs 
which appear in later chapters. ‘This, though it weakens the logical structure 
of the book, is not a bad feature in what is obviously a textbook designed to 
buttress a parallel teaching programme. ‘The first five chapters are devoted 
to elementary fluid dynamics, with illustrative examples. ‘lwo more 
specialized chapters on pipe flow and on pumps and compressors follow. 
For some reason turbines are discussed in an earlier chapter, amongst the 
miscellaneous applications of the basic ideas. 

Kay next turns to heat transfer with brief treatments of conduction and 
heat exchangers. Dimensional analysis is taken up again and applied to heat 
transfer, and the co-existence of skin friction and heat transfer is discussed 
in another short chapter. ‘The subsequent six chapters are rather strangely 
described in the preface as forming a hard core of basic theory. ‘The signifi- 
cant word here presumably is ‘ hard’, since more mathematical demands 
are now made on the student, with the development of the Navier-Stokes 
equations and simple applications thereof. ‘There is a chapter on boundary 
layers, including the momentum equation, and another containing the usual 
superficial treatment of turbulent pipe flow. Successive chapters then turn 
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to mass transfer, energy equations and again to forced convection. ‘The 
book ends with five increasingly specialized chapters on compressible flow, 


open-channel flow, solid particles in suspension, packed beds and heat 


transfer with phase change. 

or mechanical engineer, the book has to discuss cursorily a large number of 
topics. In consequence, some, particularly the later, chapters become a 
motley miscellany. ‘lhe primary interest is in internal rather than external 
aerodynamics. Potential flow is mentioned only briefly, but a dispropor- 
tionate amount of space is given to well-stated external boundary layer theory, 
both in the text and an appendix. 

[In comparison with most heat transfer textbooks, the book is refreshingly 
fundamental in outlook and refrains from too many unjustified assertions. 
The tenth chapter, in which dimensional analysis is applied to heat transfer, 
is particularly fair in developing this subject, notoriously suspect from an 
undergraduate viewpoint. ‘lhe need for five basic dimensions in describing 
phenomena which do not involve the conversion of significant amounts of 
work or kinetic energy into heat and internal energy is well explained. On 
the other hand, to take an adverse example, the usual assertion that the 
Froude number is important when there is a free surface is not supported by 
clear and convincing arguments. Even in the last three highly specialized 
chapters the basic theoretical approach is preserved. In fact there is too 
little descriptive comment on the many formulae, and as a result justice is 
hardiy done to the more empirical topics like heat transfer with boiling. 
Just occasionally the author turns to extensive practical details, notably at 
the beginning of the chapter on pumps and compressors. 

Author and publishers are to be congratulated on producing a first edition 
so free from mistakes. Among the few errors 1s the unnecessary requirement 
that incompressible flow must be steady if the divergence of velocity is to 
vanish. Some reference to the idea of thermal ratio or effectiveness in heat 
exchangers seems desirable, and a clear definition of bulk or mixed tempera- 
ture is necessary. ‘here are many details which could be singled out as 
praiseworthy, such as the use of specific weight instead of density to make the 
hydraulicist’s form of Bernoulli’s equation dimensionally consistent. ‘he 
collection of advanced examples which ends the book will be useful to 
university teachers. It is interesting to find Hartmann’s analysis of the 
simple electromagnetic pump given as one of these examples. 

Summing up, one may describe Kay’s book as a useful textbook for 
general engineering undergraduate courses. For subsequent reference, 
engineers would do better to consult more advanced books with a more 
logical structure. Students of applied mathematics would also probably 
tind such books more to their liking. Nevertheless, Kay’s book has for its 
class an unusually fundamental outlook and, as remarked earlier, marks a 
beginning in the good work of persuading engineers to master vector analysis 
and other slightly advanced mathematical techniques. 


J. A. SHERCLIFI 
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